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Abstract 


A method was developed for the geometrically nonlinear analysis of the static response 
of thin-walled stiffened composite structures loaded in uniaxial or biaxial compression. The 
method is applicable to arbitrary prismatic configurations composed of linked plate strips, 
such as stiffened panels and thin-walled columns. The longitudinal ends of the structure are 
assumed to be simply supported, and geometric shape imperfections can be modelled. The 
method can predict the nonlinear phenomena of postbuckiing strength and imperfection sen- 
sitivity which are exhibited by some buckling-dominated structures. The method is 
computer-based and is semi-analytic in nature, making it computationally economical in 
comparison to finite element methods. 

The method uses a perturbation approach based on the use of a series of buckling mode 
shapes to represent displacement contributions associated with nonlinear response. Dis- 
placement contributions which are of second order in the modal amplitudes are incorporated 
in addition to the buckling mode shapes. The principle of virtual work is applied using a finite 
basis of buckling modes, and terms through the third order in the modal amplitudes are re- 
tained. A set of cubic nonlinear algebraic equations are obtained, from which approximate 
equilibrium solutions are determined. Buckling mode shapes for the general class of struc- 
ture are obtained using the VIPASA analysis code within the PASCO stiffened-panel design 
code. Thus, subject to some additional restrictions in loading and plate anisotropy, structures 
which can be modelled with respect to buckling behavior by VIPASA can be analyzed with 
respect to nonlinear response using the new method. 

Results obtained using the method are compared with both experimental and analytical 
results in the literature. The configurations investigated include several different unstiffened 
and blade-stiffened panel configurations, featuring both homogeneous, Isotropic materials and 
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laminated composite material. Results for the local-postbuckling response of stiffened and 
unstiffened panels agree well with results in the literature for moderate postbuckling load 
levels. In flat blade-stiffened panels which exhibit significant interaction of the local and Euler 
buckling modes, the method is successful in predicting the consequent imperfection sensitiv- 
ity, but the method loses accuracy as imperfection amplitudes are increased. 
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1.0 Introduction 


Thin-walled compression members are found in a large variety of structural applications. 
These include stiffened panels found in, for example, the skins and ribs of aircraft wings, the 
bulkheads of ships, and thin-walled open- and closed-section columns used in buildings and 
bridges. Most existing structural components of this type are fabricated of metallic materials; 
however the same need for strength, stiffness, and reduced weight which makes thin-walled 
stiffened structures attractive has stimulated interest in the potential benefits of using fiber- 
reinforced composite materials. 

The specific class of thin-walled compression members considered here are 
compressively loaded prismatic configurations, some examples of which are shown in 
Figure 1. These are structures (or structural components) which can be modelled as 
prismatic assemblages of finite-length plate strips, where each strip is linked at one or both 
longitudinal edges to other plate strips and is subjected primarily to in-plane loading as a re- 
sult of the global loads imposed on the structure. These configurations are generally mod- 
elled with regard to classical buckling response, where the lowest eigenvalue provided by a 
buckling analysis is assumed to represent the maximum load sustainable by the structure, 
provided that material failure has not occurred at a lower load level. 

There are often geometrically nonlinear effects which cause an actual structure to pos- 
sess strength significantly greater than or less than the strength predicted based on a buckling 
analysis. When the primary buckling mode is a "local* mode (where the buckling halfwave 
length in the longitudinal direction is of the same order of magnitude as a representative width 
between plate junctures in the stiffened structure, and the junctures between plate compo- 
nents rotate but do not translate), then often the structure can support considerable additional 
load beyond the buckling load. Such a structure is said to have postbuckling strength. When 
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Figure 1, Examples of compressively loaded prismatic linked-piate structures. 
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the primary buckling mode is the "global" (or "Euler") mode (characterized by the deflection 
of the structure in a single halfwave over the length), the structure may have some 
postbuckling strength, or may be neutrally stable. 

When the global and local buckling modes have associated theoretical buckling 
eigenvalues which are equal, the mathematically idealized structure is highly unstable in 
postbuckling. This is due to the phenomenon called "modal interaction;" the overall bending 
stiffness which resists global postbuckling motion is reduced by local postbucking defor- 
mations. When an unstable postbuckling path exists for an idealized structure, a physical 
structure generally exhibits elastic limit-load behavior, where the limit load is less than the 
theoretical buckling load. This is because imperfections and secondary loads which are inev- 
itably present (material and geometric Imperfections, load eccentricity, lateral pressure, 
structural weight, and inertial forces) make the structural response nonlinear from the onset 
of load application, causing the load-response path to shift away from the theoretical 
prebuckling path toward the unstable postbuckling path without achieving the theoretical 
buckling load. Structures exhibiting this type of modal interaction are thus said to be 
imperfection sensitive with regard to elastic strength. Because a structure may exhibit this 
type of imperfection sensitivity, or may exhibit postbuckling strength, the accurate prediction 
of the strength requires an analysis which includes both the effects of geometric nonlinear- 
ities, and the effects of imperfections. 

A great deal of emphasis is placed on mass minimization of stiffened structures, partly 
for the purpose of minimizing material usage, but more importantly for transportational vehi- 
cles, to minimize fuel consumption with a concern for life-cycle costs and environmental im- 
pact. In aeronautical and space applications, weight minimization receives great emphasis, 
so that specified safety factors for strength are often relatively small. Based on the failure 
criteria used during the optimization process, a minimum-mass design may satisfy safety 
factors without any additional margins, and thus, inaccuracies in the prediction of failure can 
result in unacceptable strength or performance characteristics, or conversely, in undesirable 
weight penalties. An optimal stiffened-panel design which is based on buckling constraints is 
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particularly suspect with regard to adequacy of strength, because such a design is prone to 
have simultaneous buckling eigenvalues for the local and global buckling modes [1 J; such a 
design is often termed a "naive" optimum design, because the possibility exists that a signif- 
icantly different design is truly optimal if the effects of imperfections are taken into account. 
Thus, it is desirable to replace or augment the classical buckling analysis with a correspond- 
ing geometrically nonlinear analysis. 

Finite element analysis techniques exist which can compute the elastic behavior of the 
described structures to practically any desired degree of accuracy, but these techniques are 
computationally expensive, particularly when nonlinear effects are modelled, and thus are 
poorly suited for applications in optimization or preliminary design. Several stiffened-panel 
design codes have been developed which use analytical or semi-analytical approaches to 
obtain more economical methods of geometrically nonlinear analysis. These include PASCO 
{Panel Analysis and Sizing Code) [2], POSTOP (Postbuckled Open-Stiffener Optimum Panels) 
[3], and PANDA2 [4], each of which attempts to account for some or all of the nonlinear phe- 
nomena discussed above. However, each of these methods has limitations with regard to the 
configurations which can be accurately modelled. 

As the subject of this dissertation, a new method has been developed for the geomet- 
rically nonlinear analysis of the equilibrium response of the class of compressively loaded 
prismatic composite structures which can be modeled as a set of linked plate strips. The 
structures are assumed to have simple support at the longitudinal ends, and may have any 
of a variety of support conditions specified at the side boundaries. The component plate strips 
may have anisotropic material properties suitable for modelling a variety of composite lami- 
nates of practical interest. Loading of the structure may be uniaxial or, for some panel-like 
configurations, biaxial, and the method can be used with or without geometric shape 
imperfections. The method employs buckling eigenfunctions, obtained from the VIPASA 
buckling and vibration analysis program [5] (which is Incorporated in the PASCO design code), 
as shape functions for modelling geometrically nonlinear response; thus, subject to some 
additional restrictions in loading and plate anisotropy, structures which can be modelled with 
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regard to buckling response by VIPASA can be analyzed with regard to geometrically nonlin- 
ear response using the method presented here. The method is encoded in a FORTRAN pro- 
gram named NLPAN, which has been implemented on a digital computer 

In the following section, literature pertaining to the postbuckling and nonlinear response 
of compressively loaded plates and linked-plate type structures is reviewed. Next, the theory 
for the new method is presented in detail. Finally the method is applied to several flat stiff- 
ened and unstiffened panel configurations in order to explore certain options in the analytical 
approach and to assess the performance of the method. 
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2.0 Review of Literature 


The review of literature here will be limited to a consideration of flat stiffened and un- 
stiffened panels to which only axial or biaxial in-plane loading is applied. Structures with un- 
symmetric laminates will not be considered. While curved panels and shells can be modelled 
approximately as assemblages of plate strips, the literature for such configurations will not 
be reviewed. The focus on the literature of flat panels does not reflect a limitation in the 
method of analysis presented here, but rather limits the scope of the literature review to 
problems of the type investigated in the verification test cases discussed in Chapter 4. The 
review is additionally limited to a consideration of structural response associated with elastic 
material behavior. 

There are three survey and summary papers in the literature which provide a good 
overview of issues relating to the response and analysis of compressively loaded flat panels. 
The buckling and postbuckling behavior of composite plates, and to a lesser degree, stiffened 
composite panels, is discussed in [6]. The buckling and postbuckling of stiffened panels (and 
shells) are discussed in [7] and [1], with an emphasis on the prediction of imperfection sen- 
sitivity associated with modal interaction. 

Some literature on postbuckled plates is first considered, because it helps to provide a 
perspective on the response and analysis of postbuckled stiffened panels and thin-walled 
columns. Observations from experimental investigations of postbuckled plates will be dis- 
cussed, followed by a consideration of some methods which have been used to analyze them. 
Next, literature will be reviewed which covers experimental investigations of the nonlinear 
behavior of compressively loaded stiffened panels and thin-walled columns, followed by a re- 
view of existing methods of analysis for those structures. 
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2.1 Postbuckled Plate Literature 


2.1.1 Experimental Observations of Postbuckled Plates 


Compressiveiy loaded rectangular plates with supported edges have postbuckling 
strength which is generally limited only by material yielding or failure. This fact has been 
known at least since the early days of aluminum airframe construction, and motivated the 
often-cited study of von Karman and associates [8] concerning the strength of postbuckled 
metallic plates. Simply supported composite plates can display significant postbuckling 
strength, as reported in, for example, [9] and [10]. Uniaxially loaded composite plates with one 
side edge simply supported and one side edge free also demonstrate postbuckling strength 
[ 11 ]. 

For a simply supported rectangular plate with orthotropic material properties, buckling 
mode shapes are sinusoidal in the directions of both of the reference axes. During progres- 
sive loading into the postbuckling regime, the postbuckled shape usually retains the same 
qualitative shape as the buckling mode. In some circumstances, however, uniaxially loaded 
plates (particularly those which are longer in the direction of the load axis than they are wide) 
exhibit a sudden increase in the number of halfwaves in the direction of the load axis [10,12]. 
The same type of mode change has been observed in one-edge-free plates [11]. In two 
nominally identical one-edge-free postbuckling specimens tested in [11], one exhibited a mode 
change, whereas the other did not. This suggests that imperfections in the plates or test fix- 
tures may play a role in producing the mode change phenomenon. The phenomenon of mode 
change in postbuckled plates appears to be benign in that it did not, in the cited experiments, 
cause detectable material failure, nor did it coincide with a limit in elastic postbuckling 
strength. However, the sudden snap observed in the experiments may cause subtle material 
damage which accumulates with repetitive snapping due to load cycling, resulting in fatigue 
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failure. Thus, the prediction of mode-change during postbuckling response may be important 
in structural applications where plates or stiffened panels are designed to operate in 
postbuckling. 


2.1.2 Analysis of Postbuckled Plates 

Levy [13] and Coan [14] provide postbuckling solutions for simply supported square 
plates which are loaded uniaxially with uniform end shortening. The form of the displace- 
ments used to represent solutions in [13] and [14] serves as a reference form, to which the 
form of the displacement solutions used in the current method (see Section 3.2.1) are com- 
pared. Thus, the solutions of [13] and [14] are examined here in some detail. The von Karman 
nonlinear plate equations are assumed to apply. In Levy's study, the unloaded side edges 
are constrained to remain straight with no net normal in-plane loads, whereas in Coan's study, 
the unloaded side edges are uniformly free of in-plane loads, and thus they may experience 
in-plane distortion. Both studies express out-of-plane deflections w as a double Fourier sine 
series (where the notation of [14] has been converted to the notation of [13]): 

OO OO 

mn sin(mn-x/L) sin(mry/b) 0) 

m = In = 1 

where w mn is a load-dependent scalar, x and y are the reference coordinates in the plane of 
the plate, L is the dimension in the x-direction (the load axis), and b is the plate dimension in 
the y-direction. Renaming the mn th load-dependent constant, w mn = <7„ and shape function 
sin(mrrx/L) sin (nny/b) = w,(x, y), the complete displacement field for the postbuckled solutions 
of both [13] and [14] can be expressed in the form 
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(2) 


u(x, y) = Au l (x) + Y, Z W/ x ' ^ 
/= ly = 1 

oo oo 

v(x, y) = ;.^(y) + Z Z W(/< x - y) 
/ = i y = i 

oo 

w(x, y) = Z<W( x ,y) 


; = i 


where u and v are the in-plane displacements in the x- and y-directions, respectively, X is the 
parameter controlling end shortening, Xul and Xvl give the displacements for the unbuckled 
solution, and u,,(x, y) and v /; (x, y) are load-independent shape functions. 

There are several points worth noting regarding the solution form shown above. First, 
each function w, represents a buckling eigenfunction for the plate, and thus, the load- 
dependent scalars g,- represent "modal amplitudes" for the participating buckling mode 
shapes. Second, when the above form is used to express the nonlinear equilibrium equations 
for the plate, each If pair of in-plane shape functions, u u and can be determined inde- 
pendently to identically satisfy the in-piane equilibrium equations and the homogeneous form 
of the boundary conditions; thus, in truncating the infinite series expression for w in order to 
make the problem tractable, the in-plane equilibrium equations are always satisfied. The third 
point noted here is that retention of only a few terms in the series for w (approximately four 
to six) provides results for displacements which are fairly accurate up to loads of several 
times the buckling load [13]. The primary buckling mode provides the dominant contribution 
to the out-of-plane displacements in postbuckling. It is also interesting to note that while the 
difference in the in-plane boundary conditions used by Levy versus those used by Coan does 
not affect the buckling load or mode shape, there are significant differences in the 
postbuckling response (see Section 4.1). This indicates a sensitivity of postbuckling response 
to variations in the boundary conditions at the side edges. 

The literature contains many additional works dealing with the postbuckling analysis of 
plates; a large variety of plate shapes, anisotropic properties, and boundary conditions are 
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considered. Many of these works are summarized in [15]. The great variety of analyses in the 
literature are not discussed further here, except for the topic of predicting mode changes in 
postbuckling. Mode-change is investigated in [16], where out-of-plane displacements are re- 
presented in terms of only two buckling modes; it was found necessary to include the effects 
of shape imperfections in order to model the mode-change phenomenon. Direct minimization 
of the total potential energy is used in [17] to compute equilibrium solutions and predict the 
mode-change phenomenon. Mode change is assumed to be consistent with the assumption 
that the plate seeks the absolute lowest total potential energy level. This method does not 
rigorously model the response of a plate to a monotonically varied generalized load. A stable 
equilibrium state can exist at any local minimum in the total potential energy, so that the 
change in mode shape actually occurs when the active local minimum disappears with in- 
creasing load. This discrepancy is avoided in the theoretical approach of Thurston [18], which 
is capable of tracing the load-response path through limit points and secondary bifurcation 
points. 


2.2 Postbuckled Stiffened-Panel Literature 


2.2.1 Experimental Observations of Postbuckled Stiffened-Panels 


As was the case for simple plates with supported edges, compressively loaded stiffened 
panels can exhibit considerable postbuckling strength. This has been observed in, for exam- 
ple, blade-stiffened composite panels [19], l-stiffened composite panels [20], hat-stiffened 
composite panels [19,21], and in simple one-edge-free, channel, hat, I, and Z composite 
stiffener sections [11,22], The significant postbuckling strength observed for the structures 
listed above is primarily associated with buckling and postbuckling response which occurs in 
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a "local" mode shape, i.e. one in which the buckling mode shape is generally characterized 
by rotation of plate-strip junctures without translation, and by a buckling halfwave length in the 
longitudinal direction which is of the same order of magnitude as a representative width be- 
tween plate junctures in the stiffened structure. 

A set of three J-stiffened composite panels were tested for buckling response in [23], The 
panels were of a variety of lengths, but each had two stiffeners and unsupported side-edges, 
and the loaded ends of each panel were potted. All three of the panels buckled in a "torsional" 
(or "twisting") mode, characterized by rotation of the stiffener sections about the lines of at- 
tachment to the panel skin, with compatible rotation of the skin at the joint line. The modes 
occur in halfwave lengths much greater than what would be expected for a pure local buckling 
mode in the sense described previously. Two of the three panels tested did exhibit 
postbuckling strength (the third suffered material failure upon buckling), but for each panel, 
the observed buckling load was significantly less than the theoretical buckling load. It was 
hypothesized that the panels are imperfection sensitive with regard to the onset of the twisting 
buckling mode [23]. 

The global (also called Euler or wide-column) buckling/postbuckling mode is character- 
ized by a bowing type of response occurring over the length of the panel with little distortion 
of the cross-section, assuming that the side-edges are unsupported. If the side edges are 
supported, then the panel will bow in two coordinate directions in the manner of a simply- 
supported plate. The global mode of response is not investigated much in the recent literature, 
because it can be analyzed rather easily using wide-column theory or plate theory in con- 
junction with the extensional and bending stiffnesses of the panel cross-section. Some ex- 
perimental results are presented by Tulk and Walker [24] for a wide blade-stiffened panel of 
homogeneous, isotropic material with simple end supports and unsupported side-edges, 
subjected to uniform end-shortening. For cases where no local buckling is induced, the 
postbuckling response is seen to be essentially neutral, meaning that for end-shortening val- 
ues beyond the critical value, the load carried by the panel remains at about the critical value. 
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As is the case for rectangular plates, sudden changes in the mode shape during local 
postbuckling are observed for stiffened structures. A sudden increase in the number of 
halfwaves in the direction of the load axis is observed for locally postbuckled composite 
channels in [11], Two nominally identical channel specimens were tested in compression in 
[11], One specimen exhibited a mode change, whereas the other did not, again suggesting 
that imperfections in the specimens or test fixtures may play a role in producing the local- 
mode-change phenomenon. 

The phenomenon of mode-change discussed above is often referred to as a manifestation 
of "modal interaction," because the two different displacement patterns observed are similar 
to two of the buckling mode shapes. While this form of modal interaction involves a sudden 
switching between two dominant shapes, there is a second type of modal interaction exhibited 
by stiffened panels which has received considerable attention since the mid-1960's. This 
second type of modal interaction occurs between local buckling and global buckling modes 
having critical loads which are relatively close in value. This type of modal interaction can 
cause the elastic strength of the structure to be highly sensitive to geometric imperfections. 

in order to study the interaction of local and global buckling modes, a series of tests were 
performed on wide blade-stiffened panels of homogeneous isotropic material [24-27], All 
panels had simply supported ends and unsupported side-edges. Panels with two different 
cross-section types were considered. The first type had stocky stiffeners which deformed little 
during local buckling and postbuckling, and the second type had thin stiffeners which partic- 
ipated significantly in local buckling and postbuckling. Panels of a wide variety of lengths were 
tested, and the effects of imperfections in the shapes of both the local and global buckling 
modes were studied. Imperfection sensitivity was demonstrated to be severe in some cases, 
causing as much as a 44% reduction in the elastic limit load compared to the theoretical 
buckling load. The sensitivity was found to be most extreme for configurations in which the 
local and global buckling eigenvalues were approximately the same, but significant 
imperfection sensitivity was observed in panels for which the two buckling eigenvalues dif- 
fered by as much as a factor of two. 
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2.2.2 Analysis of Postbuckled Stiffened-Panels 


2.2.2.1 General comments 

The computer code STAGSC-1 [28,29] appears to be one of the most powerful tools 
available for the general nonlinear analysis of shell-type structures (including stiffened and 
unstiffened panels). It can use a variety of structural models, including finite elements, and its 
greatest strength is its implementation of strategies for dealing with complicated nonlinear 
phenomena. For example, limit points can be traversed using the constant-arc-length method 
of Riks [30], and the method of Thurston [31] is available in the code in order to navigate 
bifurcation off of a nonlinear load/response path, including the case where multiple bifurcation 
paths coincide [32], STAGSC-1 represents a superior analytical method for assessing results 
obtained using simplified nonlinear analyses, but it is extremely expensive when used to 
simulate the nonlinear behavior of complicated structural configurations. The approaches to 
analysis reviewed in this section attempt to model nonlinear phenomena by more economical 
means. While the issues of traversing limit points and bifurcation points are important in 
performing a general nonlinear analysis, the emphasis in this review of analysis techniques 
is on the formulation of the solution to the governing equations as it affects the accuracy of the 
solutions obtained. 

In 1945 Koiter [33] presented a general theory of the stability of equilibrium for elastic 
systems with conservative loading, and his work has served as the basis for a majority of work 
done since. Although Koiter addresses "stability," he addresses also the effects of geometric 
imperfections and the phenomena of postbuckling strength and imperfection sensitivity, all of 
which are demonstrated to have important practical consequences in structures for which 
stability is an issue. Koiter used a perturbation approach to study the equilibrium behavior 
of perfect and imperfect structures in the vicinity of the theoretical bifurcation point. While the 
approach of [33] has the strength of facilitating the analytical determination of parameters 
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which characterize the initial postbuckling response, the method is limited to a consideration 
of structural response in the neighborhood of the critical bifurcation point. Nonetheless the 
following important conclusions are drawn from [33]: i) buckling eigenfunctions can serve as 
useful shape functions in the nonlinear analysis of buckling-dominated structures, by helping 
to minimize the number of independent variables in an analysis, ii) the phenomenon of 
postbuckling strength is associated with a stable secondary equilibrium paths, iii) the phe- 
nomenon of Imperfection sensitivity is associated with an unstable secondary equilibrium 
path, and iv) the degradation of load-carrying capability in imperfection-sensitive structures 
is dependent on the amplitude, and sometimes the sense (positive or negative), of shape 
imperfections. Regarding the third point above, it has been demonstrated analytically that for 
blade-stiffened panels with coincident critical loads for local and global buckling, a highly un- 
stable postbuckling equilibrium path exists [34], This confirms a theoretical basis for the 
possibility of imperfection sensitivity in stiffened panels. 

The greater role played by transverse (out-of-plane) shear deformation in the buckling 
of composite plates compared to metal plates is well established [35], and thus, some authors 
include transverse shear deformation effects in buckling and postbuckling analyses [35-37], 
Similarly, the nonlinear material behavior exhibited by some composites is emphasized by 
some authors in the analysis of composite structures [36]. Nonetheless, both transverse shear 
deformation and nonlinear material behavior produce higher-order effects which it seems 
reasonable to ignore for the sake of simplicity in developing a new method such as the one 
presented here. 

One aspect shared by nearly all analytic or semi-analytic (mixed analytic and numerical) 
methods of analysis for postbuckling is the assumption that buckling mode shapes have har- 
monic forms along the length of the component plates. This assumption reduces the partial 
differential equations governing buckling to ordinary differential equations in terms of the 
transverse in-plane coordinate [5]. Similar simplifications can be made at many other steps 
in the mathematical derivations which lead to a method of nonlinear analysis, as will be shown 
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in Chapter 3. The sacrifice in obtaining these simplifications is a restriction in the type of 
support conditions which can be simulated at the longitudinal ends of a panel. 

Based on the different elastic response modes observed in experiments, the analyses 
found in the literature can be discussed in terms of three distinct classes of behavior: global 
postbuckling, local postbuckling, and modal interaction. These classes are each considered 
in following sections. 

2.2.22 Analysis of global postbuckling 

In the analysis of global buckling and postbuckling, many researchers use the Bernoulli 
kinematic assumptions in performing wide-column analyses of stiffened panels [3,4,38], The 
Bernoulli assumptions for thin-walled prismatic structures are analogous to the Kirchhoff as- 
sumptions for beams. According to the Bernoulli assumptions, cross-sections of a member 
remain flat and perpendicular to a longitudinal reference axis which deforms with the member 
in global buckling. These kinematic assumption, while valid for some configurations, will be 
unconservative for others because in-plane shearing of stiffeners and webs is not accounted 
for, and are not applicable to panels which are of finite width with supported side edges, for 
which the cross-sections must deform during global buckling deflection. 

Giles and Anderson [39] established beam-column relationships for application to wide 
columns, which relate the mid-length bending moment to the axial load, pressure loading, in- 
itial bowing imperfection, and axial load eccentricity. In the PANDA2 design code [4], two 
options for a global buckling analysis are available. The first option is to represent a finite- 
width stiffened panel as a finite-width unstiffened panel with smeared stiffener properties, al- 
lowing the use of plate theory. The second option is a wide-column buckling analysis for a 
representative stiffener/skin unit of the panel. Bushnell [4] and Anderson and Stroud [40] 
simulate wide, clamped-end panels by modelling reduced-length panels with simply supported 
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ends. Bushnell uses a panel length reduced by the divisor ,/3.85 ; Anderson and Stroud select 
a divisor value of 2. 

22.2.3 Analysis of local postbuckling 

Local buckling analyses often use methods along the lines of plate buckling analyses, 
with the added complication that displacement compatibility and equilibrium must be satisfied 
where plates join. This has the consequence that the load multiplier is no longer separable 
as an isolated factor in the eigenvalue problem [5J, so that a search technique must be used 
to locate the buckling eigenvalues. Assuming that displacements in a nonlinear analysis are 
represented as a perturbation expansion in the buckling modal amplitudes, the higher-order 
displacement contributions associated with nonlinear response are also more difficult to de- 
termine when compared to simple plate problems. For simple plates, the displacement con- 
tributions which are of second order in the modal amplitudes involve only in-plane 
displacements (see equations (2)), but for linked-plate problems, the in-plane displacements 
of one plate strip induce out-of-plane displacements in (noncoplanar) adjoining plate strips. 

To reduce these complications, some authors use what Sridharan and Peng [41] refer to 
as the "classical assumptions" of local postbuckling deformations. These are [41]: "1) at a 
corner of the plate structure where two plates meet at an angle, the normal displacement w 
for each plate vanishes; 2) the normal stress resultant N y in the transverse direction vanishes 
for each plate at the corner." A quantitative justification for the assumptions was given in [42]. 
For a case such as a blade-stiffened panel where one plate (the stiffener) joins at an angle to 
a second plate which is locally continuous (the panel skin), the "classical assumption as 
stated can be modified to improve the physical accuracy of the assumptions. In the modified 
form, the stiffener introduces no discontinuity in the normal stress resultant N y in the panel 
skin. This allows for a nonzero value of N y in the skin at the line of attachment of the stiffener. 
In certain circumstances the "classical assumptions" discussed above are fundamentally in- 
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correct or inaccurate, such as when three plate strips with a common joint lie in three different 
planes, or when two adjoining plate strips lie in planes which are nearly, but not exactly, 
parallel. 


2.2.2.4 Analysis of modal interaction 

The analysis of modal interaction and imperfection sensitivity usually involves the use of 
modal amplitudes as generalized coordinates. Imperfection shapes are also expressed in 
terms of modal imperfection amplitudes corresponding to the buckling mode-shapes included 
in the analysis. In other words, the modal imperfection amplitudes are the Fourier coefficients 
of a series expression for the imperfection shape in terms of buckling eigenfunctions. 

The first study directed at analyzing the imperfection sensitivity of compressively loaded 
linked-plate structures is generally attributed to van der Neut [43]. He modelled an idealized 
thin-walled column with a rectangular cross-section. The two flanges are assumed to carry 
the entire axial load, and they are held in relative position, with simple support at the side 
edges, by two webs which are infinitely stiff in shear but which support no axial load. A 
local/global division is used in the analysis, in which the global buckling/postbuckling behav- 
ior is based on the tangent extensional and bending stiffnesses of the idealized cross-section, 
which in turn depend on the curvature due to global bending. The local postbuckling analysis 
uses analytical results for the postbuckling response of long, simply supported rectangular 
plates strips, in order to compute the axial stiffness of each postbuckled flange. The axial and 
bending stiffnesses of the column cross-section are thus easily determined. Imperfections in 
the shape of local and global buckling modes are accounted for. Van der Neut's results pre- 
dict imperfection sensitivity in the idealized column, particularly when the theoretical global 
buckling load is equal to or somewhat greater than the load for local buckling. Imperfections 
in the shape of the local buckling mode lead to global buckling at loads below the theoretical 
local buckling load. In [44], van der Neut applied a similar idealized approach to the analysis 
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of hat-stiffened panels, in which the panel skin and the cap of the hat are assumed to carry 
all of the axial load. Imperfection sensitivity was found to be less severe than in the thin-walled 
column. Crawford and Hedgepath [45] use a similar approach to the analysis of a truss-core 
sandwich panel. A certain degree of imperfection sensitivity was predicted. The general ap- 
proach of a local/global division in the analysis is adopted in many later studies. 

Around the same time as van der Neut's study, Graves-Smith [46] performed a detailed 
analysis of the interaction of local and global buckling in thin-walled metallic columns of rec- 
tangular section. Plastic deformation of the material was also accounted for. The results pre- 
dicted imperfection sensitivity in the elastic regime, due to interaction of the local and global 
buckling modes. Graves-Smith computed the deflection fields which were quadratic in the 
modal amplitudes (as were shown to play a role in postbuckled plates), of which one is bi- 
linear in the local and global modal amplitudes. This latter displacement field, which was ob- 
served experimentally by Bijlaard and Fisher in 1953 [41], has the same approximate 
cross-sectional shape as one of the secondary (higher in eigenvalue) local buckling modes 
with the same number of halfwaves along the length of the column as the primary local 
buckling mode. This observation has led to approaches in which the secondary local buckling 
mode is included in the analysis as another interacting buckling mode [38,47-49], This ap- 
proach offers certain advantages which are discussed in later sections. 

A blade-stiffened panel configuration studied by Tvergaard [50,51] has subsequently been 
studied by several other researchers [4,38,52-54]. The Tvergaard panel is a wide panel with 
uniformly spaced blade stiffeners (without flanges) on one side of the panel. The panel is 
subjected to uniform end-shortening at simply supported ends (loaded along the neutral 
bending axis). The configuration is classified as singly symmetric in the two dominant buckling 
modes. This means that the local buckling mode is symmetric with regard to positive or neg- 
ative sense, whereas the global mode is not symmetric; global buckling deflections of one 
sense will increase the compression on the skin and decrease the compression near the 
stiffener tips, whereas global deflections of the other sense will have the opposite effect. In 
reference [50], Tvergaard used Koiter's perturbation approach of [33] to analyze a configura- 
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tion with coincidental buckling eigenvalues for the local and global buckling modes. The con- 
tributions to deflection fields were limited to the linear in-plane displacements, and 
out-of-plane displacements in the shapes of the two buckling modes. Total potential energy 
terms through the third order were retained. An often-cited result of Tvergaard's (Figure 4 of 
[50]) is the fall-off in the limit load with increasing imperfection amplitudes in the shapes of the 
local and global buckling modes. For the test case reported, imperfection amplitudes which 
are on the order of the plate thickness cause a reduction in load capacity of over fifty percent. 

Subsequent, more accurate analyses by Tvergaard [51], Koiter and Pignataro [52], and 
others have led to the conclusion that while the results in [50] are asymptotically correct as 
the imperfection amplitudes go to zero, the results are overly pessimistic for imperfection 
amplitudes of practical significance. The analyses of [51] and [52] both feature the following 
significant differences from the original analysis of Tvergaard: potential energy terms are re- 
tained through the fourth order in the modal amplitudes rather than the third order, and dis- 
placement contributions which are quadratic in the modal amplitudes are included. This 
implies that for solutions of reasonable accuracy, it may be necessary for one or both of the 
above-mentioned elements to be present in the analysis. 

Amplitude modulation: Koiter and associates [47,52,53,55] used a concept called "ampli- 
tude modulation" in order to provide simple, relatively accurate analyses for thin-walled col- 
umns, and wide panels with evenly spaced stiffeners. Its use was inspired by the fact that the 
interaction of the local buckling mode (which features multiple halfwaves of uniform amplitude 
along the length) and the global buckling mode (which introduces a bending curvature that 
varies as a simple half sine-wave over the length) cause the local buckling waveform to be 
modulated, meaning that the amplitudes of the wave peaks are caused to vary along the 
length of the member. Similarly, the waveforms of all additional displacements induced by the 
interaction of the two mentioned modes exhibit amplitude modulation. The representation of 
an amplitude-modulated local buckling mode using unaltered buckling modes requires se- 
veral modes of similar transverse profile, but with different halfwave lengths. 
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In order to model the phenomenon of amplitude modulation without resorting to the use 
of a large number of buckling modes, Koiter and Kuiken [55] introduced the "amplitude mod- 
ulation function," which is used as a multiplying factor for local buckling mode shapes. To il- 
lustrate this concept, assume that x is the longitudinal coordinate axis, and that the 
x-dependence of a local buckling mode is given by sin(/THrx/L) where m is the integer halfwave 
number and L in the length of the structure. Using the modulation function f(x), the modulated 
local buckling waveform will have an x-dependence given by f(x) sin(m«x/L). If more than one 
local buckling mode is used as a shape function, then each mode has its own amplitude 
modulation function. The modulating functions are initially unknown and are determined as a 
part of the solution procedure. A benefit to the use of the amplitude modulation concept is that 
each of the participating modes represents a distinct type of structural deformation. This 
contributes to an intuitive understanding of the structural response. A drawback of the am- 
plitude modulation concept as used in the noted references is that it involves the assumption 
that the global buckling response is one-dimensional, as is the case for columns or wide col- 
umns (wide panels with unsupported side edges). It is not clear whether the concept can be 
meaningfully applied to configurations for which the global buckling deformations are highly 
two-dimensional, such as would be the case for a panel with supported side-edges. 

While the early studies of Koiter et al which employed amplitude modulation concepts 
[52,53,55] made use of only the global buckling mode and the primary local buckling mode, it 
is shown in [47] that when significant local buckling deflections occur on both sides of the 
neutral axis, the method must include a second local buckling mode. The second local 
buckling mode is one which has the same longitudinal halfwave number as the primary local 
buckling mode, but which has a transverse profile similar to that of the so-called mixed 
second-order field associated with the interaction of the global mode and the primary local 
mode. The second local buckling mode is represented in conjunction with an amplitude 
modulation function. 

Doubly symmetric structures: In the earlier discussion of the Tvergaard panel, it was 
noted that the structure is singly symmetric. It is shown in [56] that for structures which are 
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doubly symmetric with regard to buckling (such as a thin-walled rectangular column, or a 
panel stiffened symmetrically on both sides), the omission of total potential energy terms 
which are of order greater than three in the modal amplitudes {as is done in many applications 
of Koiter's general theory) precludes the prediction of modal interaction. Thus, a general 
method of analysis must consider the equivalent of total potential energy terms through the 
fourth order. 

2.2. 2. 5 General methods of analysis 

Probably the most highly developed methods for the study of imperfection sensitivity and 
modal interaction in compressively loaded linked-plate structures are those of Sridharan and 
his associates, late versions of which are described in [38,49]. The methods make use of many 
approximation concepts developed by Koiter et al [47,52,53,55]. The methods are restricted 
to columns of general cross-section or to wide panels with evenly spaced stiffeners, with ho- 
mogeneous, isotropic material properties. The Bernoulli kinematic assumptions {Section 
2.2.2. 2) are used to describe global buckling mode-shapes, and amplitude modulation func- 
tions {Section 2. 2.2.4) are employed. The methods thus require the use of only three or four 
buckling eigenfunctions to obtain accurate predictions of elastic limit loads for imperfect 
structures featuring modal interaction. A finite strip method is used for determining local 
buckling mode-shapes and second-order displacement contributions, and a finite element 
method is used to discretize the amplitude modulation functions along the longitudinal axis 
of the structure. The "classical assumptions' of local postbuckling analysis, discussed in 
Section 2. 2. 2. 3, are used in appropriate steps of the analysis. A drawback of this general 
method is the limitation imposed by the Bernoulli assumptions, namely that the global 
buckling response is one-dimensional. 

The PASCO design code for stiffened panels [2] has the capability to account for some 
geometrically nonlinear effects. The heart of the PASCO design code is the VIPASA computer 
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code [5] for the analysis of the buckling and vibration of stiffened composite panels, so the 
method of VIPASA merits attention here. VIPASA computes buckling (and vibration) 
eigensolutions for prismatic linked-plate structures subjected to in-plane loads, using a finite 
strip method. The displacement functions used on each plate strip satisfy the buckling 
equations exactly, and along the lines where plate strips join, the conditions of equilibrium 
and displacement compatibility are satisfied. Thus, the solutions are considered to be exact 
within the context of non-shear-deformable plate theory, and within the limitations of the 
boundary conditions at the panel ends. An important feature of VIPASA is its ability to obtain, 
with certainty, any specified number of eigensolutions at a given longitudinal halfwave-length, 
beginning with the one corresponding to the critical eigenvalue, and continuing with solutions 
corresponding to the ascending sequence of eigenvalues. 

The panel sizing code PASCO [2] which uses the VIPASA analysis makes use of the. 
beam-column relationships derived in [39] to determine a uniform bending moment which is 
applied to a panel to account for the effects of bowing imperfections, pressure loading, and 
eccentricity of the end load. The uniform bending moment has the effect of redistributing the 
compressive stresses in the panel section, so that the local buckling mode-shape and the lo- 
cal buckling load will be modified. The global buckling analysis is not affected. While this 
simple correction is a first step toward accounting for imperfections, the effects of 
imperfections in the shapes of local buckling modes are not modelled. An additional limitation 
is that the beam-column relationships apply only for configurations for which the global 
buckling response is one-dimensional. 

The two panel sizing codes POSTOP [3] and PANDA2 [4] use analyses with some simi- 
larities of approach. Both use a local/global split of the analysis. A panel is assumed to have 
many evenly spaced stiffeners, thus permitting use of a representative unit-width panel cell 
when analyzing local buckling and postbuckling. The local postbuckling deformation is de- 
pendent on the value of curvature at the mid-length due to global postbuckling (overall bow- 
ing). Deformation in the shape of the global buckling mode is based on the beam-column 
relationships mentioned earlier and depends on the tangent extensional and bending stiffness 
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properties of the panel, as determined by the local postbuckling state. Thus, the solution 
process iterates between the local and global levels. Both codes are limited to panels which 
consist of a continuous skin with attached stiffeners (or, in the case of PANDA2, truss-core 
sandwich panels [57]), with POSTOP additionally limited to configurations with open-section 
stiffeners. The two design codes use different approaches to the analysis of the unit cell, but 
these approaches are not discussed here. There are other significant differences in the 
modelling capabilities of the two codes. POSTOP can model only flat stiffened panels, and can 
not model the effects of an imperfection in the shape of the local buckling mode. PANDA2 can 
model cylindrical stiffened panels, and can treat imperfections in the shape of the local 
buckling mode. Because they use the unit cell concept, both POSTOP and PANDA2 are unable 
to consider the effects of support at the side boundaries, although in the absence of modal 
interaction PANDA2 can perform a smeared-stiffener analysis of a panel with supported side 
boundaries. PANDA2 makes available a large number of simple models for analyzing stiff- 
ened panels, giving the analyst many options, but requiring careful judgement as to which, if 
any, of the models are valid for computing the response of a given configuration. 

A procedure for the postbuckling analysis of laminated composite stiffened panels is 
presented by Sheinman and Frostig [58]. Blade, T, and L shaped stiffeners are accommodated 
in the procedure. The method is mixed analytic and numerical, using beam eigenfunctions to 
describe the longitudinal variation of out-of-plane displacements, and using a finite-difference 
representation to describe the distribution of the beam eigenfunctions on the cross-section of 
the panel. It is difficult to make a general assessment of the method based on the limited 
results presented in [58], 

Analytically determined buckling mode-shapes for a stiffened panel with a rather com- 
plicated cross-sectional shape are presented in Figure 1 of [59], Some of the mode shapes, 
obtained at a variety of longitudinal halfwave lengths, are unusual and defy obvious classi- 
fication as local, global or twisting modes. Some of the general methods of analysis that were 
discussed in this section rely on the classification of buckling mode-shapes Into the specified 
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categories, and thus can not be used to predict certain types of behavior which are exhibited 
by structures with complicated cross-sections. 


2.2.3 Current Research Effort 

A number of general methods of analysis have been developed which can analyze in 
economical fashion (compared to finite element analysis) the nonlinear response of 
compressively loaded prismatic structures (Section 2. 2. 2. 5). Each of these methods places 
various limitations on the range of structural configurations and material anisotropy that can 
be modelled, and two of the methods lack the ability to model local-mode shape imperfections. 
The new method described in this dissertation offers analysis capabilities which duplicate 
many capabilities of existing methods, but the new method should also add new capabilities, 
and perhaps offer a more economical alternative to some of the existing methods. 

In the current method, it is assumed that nonlinear plate theory governs the behavior of 
each plate strip in the linked-plate structure. At the heart of the theoretical approach used 
here is the expression of postbuckling displacements as an expansion in terms of buckling 
eigenfunctions. This general approach has similarities to both the series solution approach 
of Levy [13] for postbuckled plates, and the general approach to stability analysis of Koiter 
[33], Unlike the plate problem considered by Levy, the class of structures considered here 
do not, in general, admit exact series solutions in postbuckling, so a perturbation approach is 
suggested. Unlike the perturbation approach of Koiter, it is desired here to allow for the in- 
corporation of an arbitrary number of buckling mode shapes having associated critical 
buckling loads not necessarily in the proximity of the primary buckling load. The perturbation 
approach used here has the following features: 

1. The set of perturbation parameters is the set of 'modal amplitudes" corresponding to a 

set of buckling mode shapes. 
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2. In the perturbation expansion of displacements, contributions through the second order 
in the modal amplitudes are retained. 

3. A perturbation expansion of the plate Euler equations is used to obtain both the equations 
governing the buckling eigensolutions and the equations governing the second-order 
displacement fields. 

4. Any selected number of buckling mode shapes may be incorporated in an analysis. 

5. The virtual work statement for the structure is used to form a set of nonlinear algebraic 
equations governing approximate equilibrium solutions; terms are retained in the 
equations consistent with total potential energy contributions through the fourth order in 
the modal amplitudes. 

The importance of incorporating second-order displacement contributions was discussed 
in Sections 2.1.2 and 2.2.2 4, and the importance of fourth-order total potential energy con- 
tributions was discussed in Section 2. 2. 2.4. Indeed, the features described above are sufficient 
to obtain exact series solutions for the postbuckling behavior of simple rectangular plates, 
such as the solutions given in [13,14], While the current method is technically a perturbation 
method, it is hoped that with the features described, the current method will be able to provide 
accurate solutions well into the postbuckled regime, for arbitrary configurations from the 
general class under consideration. 

A major foundation for implementing the method outlined above, and a major inspiration 
for the overall effort documented here, is the VIPASA computer code [5], which is capable of 
determining the buckling eigensolutions for the general class of structures under consider- 
ation. The method of VIPASA can account for the elastic properties of some classes of lami- 
nated composite plates, and has great flexibility with regard to the configurations which can 
be modelled. Some aspects of the VIPASA analysis are discussed in later sections. A goal 
here was to develop a nonlinear analysis capability which retains the flexibility of the VIPASA 
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method; however, some additional restrictions were adopted here beyond those of VIPASA. 
While VIPASA is capable of modelling both panels with applied shear loading and panels 
whose component plates feature bending-twisting coupling, these two capabilities are not in- 
cluded in the current method. When modelled with VIPASA, both shear loading and bending- 
twisting coupling cause the buckling nodal lines to be skewed relative to the specified 
boundaries at the longitudinal ends of a rectangular panel. This brings into question the va- 
lidity of the buckling results. The computer code VICON (Vlpasa with CONstraints) [60] has 
been developed to rectify the shortcomings of VIPASA with regard to the skewing of nodal 
lines, but the buckling elgensolutions take on a much more complicated mathematical form 
than those used in the current effort. It was judged to be prudent to adopt the restrictions 
mentioned above for a first attempt at developing a general nonlinear analysis capability. The 
family of configurations which can be modelled by the current method is described in detail 
in the next chapter. 
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3.0 Method of Analysis 


This chapter begins with a specification of the defining characteristics of the stiffened 
structures to be considered, and a development of the equations which govern the equilibrium 
of the structure. The method of solution is a perturbation approach. A general form is as- 
sumed for the expression of displacements, and three different boundary-value problems are 
established which govern the shape functions used in the displacement expressions. The 
procedures used to solve the three boundary-value problem types are discussed in detail. 
The theorem of virtual work is used to obtain a set of nonlinear algebraic equations from 
which approximate equilibrium solutions are obtained. In this chapter, the term "structure" is 
used to refer to arbitrary configurations from the class of structures for which the method Is 
intended. The term "panel" is used when a particular concept under discussion, such as 
biaxial loading, is most applicable to nominally flat stiffened- or unstiffened-panel configura- 
tions. 


3.1 Equilibrium Condition and Boundary-Value Problem 


Geometrically nonlinear plate theory is assumed to govern the equilibrium behavior of 
each plate strip within a linked-plate structure, and thus this section begins with a discussion 
of the geometrically nonlinear plate theory. Next, the characteristics of a general linked-plate 
structure are discussed. Finally, equations governing the equilibrium of the structure are de- 
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veloped, both in the form of a boundary value problem, and in the form of a virtual work 


statement. 


3.1.1 Plate Strips 


The development of the equations governing the equilibrium of the individual plate strips 
is presented in detail in Appendix A. The guiding philosophy for the plate theory is presented 
in this section, and the governing equations are summarized. 


3.1.1 1 Strain-displacement equaUuns 

A component plate strip of a linked-plate structure is shown in Figure 2. The plate strip 
is rectangular in planform, and has an associated x-y-z coordinate system, where the x-y 
plane lies at the mid-surface of the undisplaced plate. The three displacement components 
of the mid-surface are denoted by {u} = [u(x,y) v(x,y) w(x,y)] T , corresponding to the x-, y-, 
and z-directions, respectively. The plate strip has dimensions of length L and width b in the 
in the x- and y-directions, respectively, with the origin of the coordinate axes located at one 
corner of the strip. Initial shape imperfections are represented using the convention that for 
an unloaded imperfect plate-strip the displacements degenerate to the Imperfection shape 
{u} = {u 0 } = [u° v° w°] r , and all stress resultants are zero. 

The in-plane strain components, {e} = [e*eyy,y] T , are restricted to be small compared to 
unity, but rotation amplitudes are permitted to be moderately large, including in-plane ro- 
tations. (Unlike plates with more conventional boundary conditions, a plate strip in a linked- 
plate configuration can undergo significant in-plane rotation if, for example, it has the role of 
a stiffener web which is participating in the global buckling of a stiffened panel.) The plate 
strip Is assumed to have a uniform thickness which is sufficiently small relative to the width 
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and length such that the Kirchhoff-Love kinematic assumptions accurately describe the dis- 
tribution of the in-plane strain components through the thickness of the plate. Given the above 
assumptions, the expression for the distribution of in-plane strains in a perfect plate strip is 
given by (equations (A5) ) 

(>(x, y, z)} = {e c (x, y)} + z(k c (x, y)} (3) 


where the mid-surface strains, {e c } = [«5 cf ySy] r , and the mid-surface curvatures, 
{k c } = [kJ K c y KSy] r , are given by (equations (A3) and (A6) ) 


ru.x + -5(4 + 

f W ’XX 1 


j v, y + .5(u, y + w 2 y ) V 

{,c c } = - 1 W, yy l 

(4) 

lu, y + v, x + w, x w, y J 

l2w, xy ) 



where, for example, u, x denotes the partial derivative of u with respect to x. 

In order to account for geometric shape imperfections, define mechanical strains and 
curvatures of the mid-surface, {e m } and {K m }, respectively, to be strains and curvatures which 
go to zero with the removal of all loads. They are given by the expressions 


'[if. 


+ ,5(v,^ + w,*)] - [u,° x + .5(v,°* + w° ! 


f rn-» 

{e } = 


{K m } = 


5 {u,y + w, y )] - lv°y + .5(u° + w ; 
[u, v + v, x + w, x w,y] - lu? y + v° + 


[^.y + 


V 1 } 

°] J 


(5) 
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xy 
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xy 


The mid-surface mechanical strains and curvatures are those to which the stress resultants 
are proportional. 


3.1 .1 .2 Conditions for equilibrium 

The plate material is assumed to be linearly elastic, and all loads are assumed to be 
applied at the plate edges. The transverse shear strains y„ and y ^ are neglected, consistent 
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with the Kirchhoff-Love assumptions, and the transverse normal stress a t is neglected. The 
theorem of virtual work for the plate strip is (equation (A15) of Appendix A evaluated for a 
rectangular plate) 


6W p = [ (N x 5e x + N y Sc c y + N xy 6y c xy + M x 6k x 4- M y S Ky + M xy dK xy ) dA 

- J o lfx* u + + U w - M x c 5vv, x ] l x = 0> L dy (6) 

- f [ f X SU + fydV -f f Z SW - MydW,y~\ | Q fi CfX 

= 0 

where A is the planform area of the plate strip, and where p denotes the plate strip index 

number. Stress resultants { N } = [/V* N y and { M } = [/Wx/Wy/Wxy] 7 introduced above follow 

standard definitions, provided in Appendix A (equations (A16) ). The edge force-resultants f x% 
f yt and f 2 are given along x-normal edges by 

f x = n x N x 

fy = n x (N xy + N x v, x ) (7) 

f z — n x (M X , X *f 2 M X y y y + N X W, X 4- N X yW,y) 

where n x = ± 1, and along /-normal edges by 


fx ~~ n y(^xy NyU'y) 

fy = HyNy ( 8 ) 

f z = ^y(2^xyx My'y “h ^xy^'X ^y^*y) 
where n y — ± 1. 

The theorem of virtual work is also expressed here in an alternate form, derived through 
the application of Green's theorem (see Appendix A equations (A24), (A27)-(A29) ): 
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( 9 ) 


6W p = - {[/¥„* + N x y,y + {NyU,y),y]6U + [ N x y ljf + Ny.y + (N ^ , x ) '^Sv 

+ [M X , XX + 2M X y, X y + My,yy + (W^W,^ + N X yW ,y), X + (A/^yVV,^ + Wy W , y ) , y 1 <5 W } (jA 

+ f [(/* - £)<5 u + (fy - fy)6v + (f z - f z )6w - (M x - M x )Sw, x l \ x = 0 L dy 
J o 

+ f [(/*- t)6U + (fy - fy)6v + (f Z - f z )5w - {My - My)6w,y]\ y= g 

J o 

= 0 

The above expression is used in establishing the plate Euler equations, and the boundary- 
value problem for the linked-plate structure. 

The Euler equations determined from the above equations are (Appendix A equations 
(A25) ) • 

N X , X +N X y,y+(NyU,y),y = 0 {^a) 

N xy'X + W y,y + (N x V,x),x = 0 ( 10fc ) 

M X< XX + 2/W x y, X y + My,yy + ( ^ X W ’ X ^ X yW,y), x + (N X yW, X + NyW,y),y — 0 (’IOC) 

The associated boundary conditions are given in equations (A26) of Appendix A, but they re- 
quire modification for application to linked-plate problems, an issue discussed in a later sec- 
tion. 


3.1. 1.3 Side-edge conditions 

The structural configurations under consideration are composed of plate strips linked to- 
gether along mutual side-edges, so the generalized displacements and generalized force- 
resultants along the side edges are given special notation to facilitate their usage at various 

A A 

steps in the analysis. The unit normal vector at a side edge is given by n = n y j , where 
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n„= 1 or -1. Each side edge of the strip is given an index number, e , e= 1, 2, and the as- 


sociated y-coordinate, y„ and y-normal unit vector component, n y , are given by 

r f [0 , -l] for e = 1 

[y e , n y-l “ ] [j, i] for e = 2 


( 11 ) 


The generalized displacements and generalized force-resultants along edge e are denoted 
by {u,} = [u e v e w e ip,] 7 , and {£} = [/» f y , f« mJ T , respectively, and follow the conventions 
shown in Figure 2. The edge rotation angle, is the gradient of w in the y-direction, evalu- 
ated at the appropriate edge: 


*•-0 *.y)l„ 


( 12 ) 


Force resultants f xe J ye , and f ze are forces per unit edge-length acting in the x-, y-, and z-di- 
rections, respectively, and m e is the edge bending moment per unit edge-length. The gener- 
alized force-resultants at the side edges are expressed in terms of the plate stress resultants 
as follows: 


^xe ~ n y(^xy N y U t y) 
fye ” ftyNy 

?ze = ^y(^^xy*x ^yy "F A/ X yW, x + NyW t y) 

m e -~ flyMy 


(13) 


3.1. 1.4 Plate constitutive equations 

The elastic properties of each plate strip are assumed to be those of a balanced, sym- 
metric laminated composite plate, where in addition the bending-twisting coupling stiffness 
terms Die and D 2 e are neglected. For this class of materials, the relationships between the 
stress resultants and the mid-surface mechanical strains and curvatures have the following 
form: 
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{I N } = [A]{e m } 


(14) 


{M} = [D]{K m } 


where 
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CM 
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Q 

1 — 1 
}> 

1 1 

II 

>A 1 2 ^22 ^ 

[D] = 

D12 D 22 0 
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o 

CD 

1 


0 0 


The above equations are a reduced form of the relationships given for general laminates in 
Appendix A equation (A30). configurations. The stiffness terms In the matrices above follow 
conventional definitions given by, for example, Jones [61], The matrices [A] and [D] may 
be different for different plate-strips within the linked-plate structure. 


3.1.2 Linked-Plate Structures 


In this section, the equations which govern the equilibrium and determine the side-edge 
conditions for each plate strip are used to specify a boundary value problem for the linked- 
plate structure. First, parameters which define the structural configuration are discussed. 
Equations are introduced for transforming quantities between the global coordinate system 
and the local coordinate systems of each plate. Conventions are given for the specification 
of panel loading and for the boundary conditions along the side boundaries of the structure. 
The nature of the boundary conditions along the longitudinal ends of the panel are discussed, 
and the conditions for equilibrium of the structure are established. 

3.1 .2.1 Specification of geometry 

The linked-plate structure has an associated global x-y-z coordinate system, where the 
global axes share a common x-scale with the local coordinate systems of the plate strips. If 
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the linked-plate structure is a nominally flat panel, then it is assumed that the global y-axis 
extends in the direction of the panel width, the overall dimensions are L in the x-direction, 
and B in the y-direction. A set of "node lines," corresponding to the free edges of plate strips 
and the joints between connected plate strips, run parallel to the x-axis. Node lines are so 
named because of the analogy with nodes in a finite element model, and are not to be con- 
fused with the node lines of a buckling mode shape. The generalized displacements at the 
side edge of a plate strip are rigidly linked to the generalized displacements of a node line. 
While each node line represents a boundary of sorts, it will be assumed in the present treat- 
ment that only two of the node lines in a configuration are designated as boundary node-lines 
upon which displacement boundary conditions or nonhomogeneous load boundary conditions 
are imposed. It is also assumed that the boundary node-lines coincide with the limits of the 
global y-domain of the panel. Some of these features and labeling conventions are shown for 
a sample configuration in Figure 3. 

Several vectors and matrices are used to establish the geometry of the linked-plate 

structure. Each plate strip is assigned an index number, p, p = 1,2 P , where P is the 

total number of plate strips. Each node line is assigned an Index number, n , n = 1,2 N , 

where N is the total number of node lines. A vector {/?„}, of length N, is defined such that the 
number of plate strips which join at node-line n is given by R n . Let Q be the maximum number 
of plate strips joining at any node line, i.e. let 0 be the largest element in the vector {/?„}. Then 
two matrices, [S„J and [£„ r ], each having dimensions N by Q, define the connectivity of the 
configuration, as follows: if p is the index number of the r 0 ' plate strip connecting to node line 
number n, and the index number of the edge through which the plate strip connects is e (see 
equations (11)), then S„, = p and E nr = e. For those nodes where fewer than 0 plate strips join, 
i.e. where R„< Q, the corresponding elements in the matrices [S„,] and [£„J are unused. 

As an example of the preceding notation, consider the stiffened panel shown in 
Figure 3, where the plate index numbers and the node line numbers are shown. The panel 
has three plate strips (P = 3) and four node lines ( N = 4), with up to three plate strips joining 
at any node line ( Q = 3). The vector { R „ ) and the matrices [S„,] and [£„,] are given by 
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Yg 


Figure 3. Stiffened pane! showing plate-strip numbers, node-line numbers, and labeling con- 
ventions: The local plate-strip axes shown don't reflect actual locations of the origins. 
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{*„} = 


( 16 ) 


[S J = 


100 

300 

123 

200 




100 
100 
2 1 2 
200 


(17) 


(18) 


where zeroes denote unused array positions. 


3. 1.2. 2 Generalized displacements and forces at a node line 

In general, there is an angular displacement n between the local y-z axes of a plate strip 
and the global y-z axes. In order to model some joint geometries more accurately, a small 
eccentricity is allowed between the side-edge of a plate strip and the associated node line. 
The eccentricity has measures e y and e z . The configuration measures are depicted along with 
the angle n in Figure 4(a), and the cross-section of an idealized configuration having nonzero 
eccentricities is shown as an example in Figure 4(b). 

The generalized displacements and generalized force-resultants along a node line, de- 
noted {U n } = [ U n V n W" 'F n ] r , and { F n } = [F? Fg F; M n ] T , respectively, are defined with respect 
to the global coordinate system, and follow the same conventions as the corresponding 
measures at the edges of a plate strip, {t/,} and {f e } (see Figure 2). In order to establish 
transformation rules relating {u,} with { U n }, consider a third coordinate system, (x y z), which 
is parallel to the local coordinate system of the plate strip, but which has it's origin located 
on the node line, as shown in Figure 4. The following transformation rules relating the gen- 
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b) Example of a configuration with nonzero eccentricities (from [5]). 


Figure 4. Eccentricity measure* e y and e„ and the local-to-global orientation angle /i. 
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eralized edge displacements in the three coordinate systems were established by Wittrick and 
Williams [5]: 


v e = V n cos ;u — W n sin ^ 
w e = V" sin n + W n cos n 

u e = u e + e z w e , x + v e , x 
v e = v e + e z tg 
Wg = w e - e y iP e 

4> e = te 


(19) 


( 20 ) 


As discussed in [5], the eccentricity transformation of equation (20) is not exact, but is suffi- 
ciently accurate if the eccentricity measures e y and e z are small compared to the characteristic 
wavelengths of the displacements in the x-direction. 

To fully evaluate the second of the above transformations, the x-dependence of the gen- 
eralized displacements must be known. Thus, the details of the second transformation will 
be discussed in later sections after the functional forms of the displacements have been in- 
troduced. The generalized force-resultants along the side edge of a plate strip, {f e }, can be 
transformed to statically equivalent values which act at along a node line and are defined with 
respect to the global coordinate system. The transformed generalized force-resultants are 
denoted {f n } = where n is the index number of the node line corresponding to 

edge number e of the plate strip under consideration. The relationships governing the trans- 
formation between { f e } and ( f n } will be discussed in later sections. 

For each node line, the generalized force-resultants {P 1 } are simply the sum of the con- 
tributions { f n } from all plate-strip edges which terminate at the node line. This relationship is 
expressed symbolically as 


{F n } = X {f\ 

r = 1 


n = 1.2 N 


( 21 ) 
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where p = S nr and e = E nr (S nr and E m are defined in Section 3.1.2.1), and subscript p indicates 
the index number of the plate strip to which the generalized force-resultants applies. 

3.1 .2.3 Panel loading 

The primary global load component applied to the configuration is the mean x-normal 
force per unit width, N*, defined as the total (tensile) end load divided by the reference width 
B. A second nonzero global load component, N y G , may be applied if the configuration is a 
panel configuration which features a continuous, initially flat skin. Component Nyc is the mean 
force per unit length (based on L) acting in the global y-direction normal to the boundary 
node-lines of the panel. In the analytical approach presented here, panel configurations fea- 
turing two or more skins could be treated by specifying the split of N# between the multiple 
skins; for simplicity, the present discussion treats only a single skin. In conventional plate 
theory, the in-plane normal stress resultants, N x and N y , are positive for tension; following this 
convention, N & and N# are both negative for compressive loads, as shown in Figure 3. 

For configurations to which biaxial loading is applicable (both N& and Nya are nonzero) 
two different options are considered for control of the load ratio. The first option is to maintain 
the load ratio NycfN , g at a constant value. The second option is to maintain the ratio of width- 
shortening to end-shortening, Av/Au, at a constant value. Both Nyc and Av are assumed to 
act across the two boundary node-lines. 

In discussing the specification of panel loading, a distinction is made between two re- 
gimes of response: 1) linear, unbuckled response of the structure ignoring shape 
imperfections, and 2) nonlinear response of the perfect or imperfect structure. In the first re- 
gime of response, it is assumed that only uniform in-plane normal strains and stress result- 
ants are induced in the component plate strips. (Whether or not this assumption is accurate 
must be judged by the analyst - see the discussion in Section 3.3.) Thus, all plate strips re- 
main free of out-of-plane displacements, and all displacements, strains, and stress resultants 
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are proportional to a load parameter X. In the second (nonlinear) regime of response, the load 
parameter X determines the panel end-shortening, and the global load components are no 
longer proportional to X, but must instead be computed after an equilibrium solution has been 
obtained. 

To initiate an analysis, a unit load N, G or a unit end shortening Au is specified, and if ap- 
plicable, a unit load Wyc or a unit width-shortening Av is specified. If the configuration does 
not feature a continuous, initially flat skin, then Av may still be specified, but N# is assumed 
to be zero for the prebuckling analysis. The specific values selected for the unit loads are 
unimportant; only the ratio of the two load components is important, and then only for the case 
of biaxial loading. Using the strategies to be presented in Section 3.3.1, any combination of 

A A 

unit input values can used to establish equivalent unit global loads N& and In the linear 
regime of response the global loads are related to the unit loads through the load multiplier 

X: 


N 


*G 


N 


Ya 




( 22 ) 


Equations (22) establishes the prebuckling load state used in the buckling analysis performed 
by VI PAS A. 


3. 1.2.4 Boundary conditions at the node lines 

Along each of the two designated boundary node-lines, with index numbers n , and n 2 , four 
boundary conditions are specified corresponding to the four degrees of freedom along a node 
line. Table 1 identifies two or three boundary-condition options for each of the four degrees 
of freedom. 
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Table 1. Options Available for Specifying Conditions Along the Boundary Node-Lines 


Option/ 

Component 

1 

2 

3 

U n ,F n x 

u n = w[ 

II 

o 

- 

C 

> 

II 

c: 

Fy — rtyNy a 

V'X = 0 

F n y=nyNy G 

w' 1 , 

W n = 0 

II 

o 

- 

V n ,M n 

-5 

It 

o 

/w" = 0 

- 


In Table 1, symbol n y ( = ± 1) reflects the unit normal vector direction at the side edges of the 
panel with respect to the global coordinate axes. Terms U[ and represent the node-line 
displacements associated with the linear response of the perfect structure to the unit load 
system. The term U[ is linear in x, corresponding to uniform axial strain. The term V f is in- 
dependent of x and corresponds to uniform displacement of the side edges {in the global 
y-direction), In option 3 for component V n , symbol F J is the mean value of the force resultant 
Ff over the length L. This boundary condition option Is useful in modelling symmetry condi- 
tions or in modelling a reinforced edge. 

To model proportional displacement between Av and Au , option 1 is chosen for compo- 
nent V n , whereas to model proportional loading between N# and either option 2 or option 
3 is chosen for component V n . If the latter option is desired, then the load Nyc will be main- 
tained in proportion to N& using the following condition: 

\/\^ = ,5 yA (23) 

where R is the constant of proportionality. 

Node lines which are not designated as boundary node-lines nonetheless 

represent boundaries with respect to the domains of the plate strips which comprise the 
linked-plate configuration. Each internal node line corresponds to either a joint line between 
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multiple plate strips, or the free edge of a plate strip. For the internal node-lines to be in 
equilibrium, the generalized force-resultants must vanish: 

n = 1, 2 N 

{F n } = {0} (24) 

n n 1 , /? 2 

3.1. 2.5 Boundary conditions at the panel ends 

Because buckling eigenfunctions provided by the VIPASA analysis [5] serve as the pri- 
mary shape functions in the proposed method of analysis, the assumptions of the VIPASA 
buckling analysis must be examined in order to understand the nature of the boundary support 
simulated at the longitudinal ends of the panel. The VIPASA buckling analysis computes "ex- 
act" buckling eigensolutions using a linear prebuckling solution in conjunction with assumed 
forms for the buckling eigenfunctions which permit the buckling equations to be satisfied ex- 
actly. The compromise associated with this approach is that the nature of the loading and 
support boundary conditions at the x-normal ends is predetermined by the assumed form of 
the buckling eigenfunctions. Consequently, the nonlinear analysis presented here has similar 
limitations with respect to the boundary conditions which can be modelled at the longitudinal 
ends of the panel. 

Corresponding to a unit solution for the prebuckling displacements [ui_v t 0], there are 
unit uniform displacement gradients, u L , x and v L , y , within each plate strip, where in addition, 
u L ,„ is constant throughout the structure, consistent with the requirement of uniform end- 
shortening. The buckling eigenfunctions, {u,} = [u, v, w,] 7 , have the following forms on each 
plate strip (a trivial phase shift has been applied relative to the conventions of [5]): 

r£ f (y) cos mjrx/L-j 

(u,(x,y)} = < tj ,(y) sin mnx/L V (25) 

'-<f>,iy) sin mnx/LJ 
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where m is a positive integer, and £,(y), r/,(y), and are known analytic functions. It can be 
seen that, in the most accurate sense, the boundary value problem is that of an infinite-length, 
uniformly compressed linked-plate structure, where the buckling displacements in the y-z 
plane (v„ w,) and the associated bending curvatures k„; and K yi are zero at 

x = nL , n = 0, ± 1, ± 2 Since x = 0 and x = L are included in these stations, then it can 

be considered that the eigenfunctions satisfy some form of simply supported boundary con- 
dition at the ends of a panel of finite length. 

The eigenfunction component u , has an x-dependence of cos m-nx/L , and thus is not 
generally zero at x = 0 and x = L. Consequently, cross-sections of the structure at x = 0 or L 
which initially lie in a plane may rotate or warp during buckling. This raises the question as 
to what boundary conditions are simulated for the u component of displacement at the ends 
of a panel of finite length. If the total displacement component u is expressed as the sum of 
the unbuckled solution and the buckling eigensolution (u = Au L + £,(y) cos mnx/L ), then it can 
be seen that over the wave length 2L\m any net shortening in x is due only to Xu L . It can 
similarly be argued that over the model length L, the effective end-shortening is due only to 
Xu L , despite the fact that for odd values of m the model length is not an integer multiple of the 
buckling wavelength, 2 L/m . Considering that any constant-x station carries the net end load, 
and that incipient buckling occurs at a constant applied end load with no addition of energy, 
there are effective axial displacements at x = 0 and L associated with energy input into the 
bounded segment of panel, and these effective end displacements do not change during in- 
cipient buckling. Thus, loading mechanisms can be imagined to exist at the ends of the finite 
length panel which duplicate the support provided by the infinite-length panel discussed, and 
the axial displacements of these imaginary loading mechanisms are governed by Xu L . 
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3 . 1 . 2.6 Equilibrium condition for the panel 


To form the virtual work statement for the complete panel, the virtual work is summed 
over all plate strips in the panel and set to zero. Stated symbolically, the virtual work state- 
ment is 


p 

Y W p = 0 (26) 

p = i 

where 6W P is given in equation (6). Substituting equations (6) into the above expression the 
virtual work statement can be expressed as 


SW = 



+ 


{M} T {SK C ))dA 



= 0 



(27) 


where the generalized force-resultants and generalized displacements along the node lines, 
{p} = [p p F; M n Y and {U"} = [U n V n W' , 'P' , ] T respectively, are now used to express the 
contributions associated with the side (y-normal) edges of the plate strips. The hats symbolize 
externally applied load components which are nonzero at points on the boundary where dis- 
placements are not specified. 

The virtual work statement is also presented in an alternate form derived by expressing 
equation (26) using equation (9): 


Method of Analysis 


45 



<5W 


K r 

-z - 

P = 1 L 


J* {C^X’X + ^xy,y + (WyW,y),y]«5U + CA/ xy , x + Afy.y + (W x V, x ). X l<5y 

+ [/W X>XX + 2M X y, X y + My,yy + (NgW^ + N^W.,,),* + ( NyyW^ + N y W , y ) , y ] <5 W } d A 


+ f C(f x - £)<5 u + (Ay - Ay)<5 v + (f z - f z )6w - {M x - M x )6w, x l l x = 0 L dy 
J o 

+ £f L ({F n }-{h) r ( dun }dx 


( 28 ) 


n = 1 


= 0 


This alternate form of the virtual work statement is used for reference in later sections. 

The Euler equations associated with the above form of the virtual work statement were 
given in equations (10), and they cause the area integral in the above equation to be iden- 
tically zero. The associated boundary conditions are obtained by considering the boundary 
integrals in the above equation along with the form of simple support at the panel ends which 
is suggested by the buckling eigenfunctions of equation (25). The boundary conditions estab- 
lished for the panel ends (x = 0, L) are 


h = h 

<5v = 0 

(29) 

6w = 0 
M x = M x = 0 


where the expression for f x is given in equation (7). The second and third conditions above 
allow for nonzero displacements associated with the linear, unbuckled solution. 

At the boundary node-lines (the side-edges of the panel), the following boundary condi- 
tions are established, where consistency with Table 1 has been imposed: 
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( 30 ) 


F* = F" = 0 or U n = XU? 

F n y = F y n = n y W y( . or V n =^ or -J-J F n y dx = F y " = n/7^ , V,£ = 0 
F” = F” = 0 or W n = 0 
M n = M n = 0 or x F n = 0 

where the notation is consistent with that used in Table 1. At the nonboundary node-lines, 
there are no nonzero loads { F n } imposed, so, repeating the statement of equation (24), 

n = 1, 2 N 

{F n } = {0} (31) 

n ¥= n v n 2 


Based on the boundary conditions established above, the virtual work statement of 
equation (27) can be rewritten as 


6W = 



+ {M} T {6K C })dA- f (f x Su)\ x = 0L dy 



dx 


= 0 


(32) 


where the displacements are assumed to satisfy the kinematic boundary conditions given In 
equations (29) and selected from (30). In the above equation, n, and n z are the node-line 
numbers for the two boundary node lines, having edge-normal vectors In the negative and 
positive global y- directions, respectively. 

A few comments are made here regarding the contribution to the virtual work from the 
boundary integral at the longitudinal ends of the panel. The longitudinal load Is imposed 
through controlled end shortening, in the solution procedure, the imposed end shortening is 
represented by the axial displacement term Au L (x), where k is the amplitude of the imposed 
end shortening, and u L (x) is the axial displacement distribution corresponding to the linear, 
unbuckled response to an imposed unit end shortening. The primary axial displacements at 
the longitudinal ends of the panel are associated with the axial displacement term Au*.(x), 
which does not contribute to Su. However the rotation and/or warping of the panel ends which 
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can occur during buckling and nonlinear response reflect the fact that Su is not uniformly zero 
at the panel ends. If it is assumed for example, that the end-shortening is applied through a 
knife-edge positioned along a neutral bending axis, then the boundary integral under consid- 
eration would contribute nothing. The boundary integral could serve as a point of introduction 
of the effects of elastic rotational restraint along the panel ends. In the present treatment the 
boundary integral under discussion will be considered to have no role. 


3.2 Form of the Solution 

3.2.1 Expression of Displacements 

Displacements are represented as the sum of linear, unbuckled contributions, and a 
truncated perturbation expansion in terms of 'modal amplitudes,' a modal amplitude being 
the amplitude multiplier for a buckling mode which is being used to represent displacements 
in the nonlinear regime of response. The modal amplitudes thus serve as the generalized 
coordinates associated with nonlinear response. Displacement contributions of second order 
in the modal amplitudes are retained. These contributions are deemed important (see Section 
2.2. 2.4) in order to achieve solutions of useful accuracy for a variety of geometric configura- 
tions. Displacement contributions of third and higher order are neglected. 

The displacements, then, have the assumed form 

{u} = k{u L ) + A{u£} + q l {uti + q l q j {u ij ) i,j = 1, 2, 3, ... (33) 

where summation over i and j is implied, and the notation {u} = [u v w] is used to refer to a 
set of compatible displacement Helds, one field corresponding to each plate strip, where the 
three components of each displacement Held are defined with respect to the local reference 
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frame of the strip. In the general expression of the approach, the indices / and j both vary from 
1 to oo, but in practice a finite number of terms are used. 

In equation (33), l is the load parameter controlling end-shortening, and A is a load pa- 
rameter which, along with A, determines the shortening of a panel across the width. Parameter 
A is required only in some circumstances, where it serves to control the ratio of the global 
load components, W^s/Nxg. in panels which admit biaxial loading (see Section 3.1 .2.3). Shape 
function {u<.} = [u t v L 0] T describes the linear unbuckled response to the prescribed unit loads 
Nxg and Nyc, which are distributed so as to produce uniform displacements of the ends and 
side edges. (Ultimately, the parameter A is replaced by an expression involving the variable 
parameters g, (/ = 1, 2, ... ) - see Section 3.6.) Shape function {u£) = [0 v£ 0] r , when used, 

A 

describes the linear, unbuckled response to a specified unit load N# = N&, where the end 
shortening held to zero. The generalized coordinate g, is the /*'’ modal amplitude, associated 
with the /"’ buckling eigenfunction, {u,} = [u, v, wj r . The second-order displacement fields are 
denoted as {(%} = [u/y v,j Wy] r . 

For an unloaded panel (?. = 0), the displacements degenerate to the imperfection shape 
of the panel: 


{u A = 0 } = {u o } = A o {u£} + q°{u l } + q?q°{u ij } i,j= 1,2,3, ... (34) 

where q? are the modal imperfection amplitudes. The form of equation (34) is selected to be 
consistent with the form of equation (33). 


3.2.2 Expansion of the Strains, Curvatures, and Stress Resultants 


The assumed forms for the displacements and the imperfection shape are used in the 
expressions for the mid-surface mechanical strains and curvatures of equations (5) to obtain 
the following expansions: 
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( 35 ) 


{ E m } = A{e l ) + (A - A°){s£ } + (q, - q ( 0 ){ 6/ ) + (qflj - q°q°){c,j} 

+ (Wflk - q°tf<lk){z,jk} + 


{* m } = (q i - + (q$i - q?qf){Ky) 


(36) 


where summation over repeated indices is implied, and the curvature and strain shape func- 
tions introduced above are given by 


{«t> - 


r*w| 

A-y 


{*?} = 


0 1 

* 


{<=/} = 


u rx 

v i'y 


L u hy 


C u ij'X "h -^( v i‘X V j'X "t" w i'x w j'x) 

{ e /j) = ^ V l]<y "h 5(U,,yU^,y + Wj,yWj,y) 

^ u l]<y "f" v ij<x -5( w rx w j<y "h w i’y w j-x) ■ 


(37) 


f v bx v jk’x + *hx w j*x "I 
i c ljk) = j u ry u jk'y + w by w jk<y j 
^ w bx w jky + *i>y w jk'X J 


C ^( v ij’x v khx + w i]<x w kbx) ") 
{Zijkt} = < - 5 ( u iry u khy + w tj'y w kky) 7 
L5(Wy, X W (f ,,y+Wy iy W W , X )J 



(38) 


In the above equations, terms have been grouped to permit the following symmetry to be im- 
posed: {u,y} = {Uji}, and { E y} = {e y ,}. This symmetry is ultimately used to reduce the computa- 
tional effort of determining the set of shape functions (u,y) (/,;'= 1, 2, ... ). Using the above 
expansions along with the plate constitutive equations (14), the stress resultants can be simi- 
larly expanded: 


{A/} = X{N l ) + (A- A 0 ){N£ } + (q, - q°){N,} + (q flj - q?q?){N$ 

+ to, <7 jQk - q?tfq°){ N ijk) + (<w/w - q°q° Q k^W,,^) 


(39) 


{M} = (q, - qf){M,} + (q,< 7 , - q?q?){M$ 


(40) 


where, for example, 
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m = 

M) = [4]{ K/ J 


( 41 ) 


and so forth, where [4] and [D] are the plate in-plane and bending stiffness matrices, re- 
spectively, with nonzero elements as shown in equations (15). The first two unit in-plane 
stress resultants are given by [N L ] = LN xL N yL 0] T and {Wf} = [A/& A/$. 0] r . The absence of the 
component N xy in these expressions, and the absence of y ^ in the corresponding strain ex- 
pressions of equations (37), are consequences of the biaxial load state and the absence of 
extension/shear coupling in the plate strips. 


3.2.3 Expansion of the Plate Euler Equations 


In order to obtain the equations which govern the various displacement shape functions 
that appear in equation (33), the plate equilibrium equations (10) are expressed using the ex- 
panded forms for displacements and stress resultants given in equations (33), (39), and (40). 
Before performing the expansion, a couple of preliminary items are discussed. First, as will 
be discussed in Section 3.3, because of the geometry of the configurations, the load- 
application details, and the elastic properties of the plate strips, the strain and stress-resultant 
expressions e xi L , N x l , e y l, A/ y i ., A/ft, c y t and Nfa are all constant over the domain of each plate strip, 
and the expressions y vL , A/*yi, cf L , yiyt, Nf yL are all zero over the domain of each plate strip. 
Consequently, some displacement contributions or their gradients are zero: 


Uuy = v L-x = u L = v L-x = ° 


(42) 


Second, as will be discussed in Section 3.6, the load parameter A can be expressed as a ge- 
ometric expansion in terms of the modaf amplitudes: 

A = q,A, + qtfjAjj + ... (43) 
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where the subscripted coefficients A„ Ay, etc. are constants. 

The plate equilibrium equations (10) are expanded and terms are grouped based their 
order in the modal amplitudes. Some resulting terms are identically zero, consistent with the 
discussion in the preceding paragraph, and are consequently deleted. For the case in which 
the imperfection amplitudes are zero (g? = 0, / = 1 . 2. 3, ... ) the following equations are ob- 
tained: 


9/CWx/x + N xyyy + /1 - W y 1 u /'yy] 

+ Qi< 1}[N Xij , X + N x y..,y + ( Ny.Ui,y),y + kNyUjJ.yy 

Qi^^xyyx + ^yyy + ^x L v i'xx^ 

+ qflj[N x y,,, X + Ny.,,y + i^X^i'X^'X ^ ^X^ij'XX 


+ AjNfu l ,yy] + 0(qf ) = 0 


+ AjN x Vj, xx ] + 0(q, ) — 0 


(44) 


d/C^x/xx + 2Mxy,,xy + ^yyyy + H^x L w i<xx + ^y L w i<yy)^ 

+ q,<7/[^x,y.xx + 2W xyjj'xy + M y u -yy + i^x^rx + N xy w,. y ), x + (N xy w hx 4- N y w ry ),y 
+ ^-(Nx L W lj’XX + NyWjj,yy) + Aj(N x W), xx + NyWj,yy)l + 0(q, ) = 0 

Bracketed expressions in the above expansions set the grounds for determination of the 
buckling eigensolutions and the second-order displacement fields. In the following three 
sections, details are given for the methods of determining the linear unbuckled solutions, the 
buckling eigensolutions, and the second-order displacement fields, respectively. 


3.3 Linear Solutions 


There are two different displacement shape-functions associated with linear unbuckled 
response of the perfect structure. The first, {u L } = v L 0], describes the response to the 
specified global unit loads N x G and where the longitudinal strain £„ is required to be uni- 

A 

form throughout the structure. As discussed in Section 3. 1.2.3, loading N& may be specified 
as nonzero only in structures in which it will not induce initial bending in any of the plate 
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strips. The second shape function, {u£} = [0 v£ 0], describes response to the specified global 
unit loading Nfc, where the longitudinal strain e x is maintained at zero throughout the structure. 
Nfc may be selected to be any nonzero value. Displacement contribution A{u£) of equation 
( 33 ) is used for configurations which have a continuous, initially flat skin, in cases where the 
ratio of the transverse to axial load is specified to remain constant. It is noted that the two 
shape functions { 14 } and {uf} are solutions to similar boundary-value problems which differ 
only by the specification of N& in the former and of e« in the latter. Thus only the solution for 
{ 14 } will be discussed, except where differences exist between the two solution procedures. 
Because of the nature of the structural geometry, the loading, and the in-plane 

A A 

constitutive equations for the plate strips (equation (15)), the unit load system N& and N# in- 
duces uniform stress resultants on each plate strip, expressed as 

N r = constant 

X L 

Ny L = constant (45) 

Nxy t = 0 

The above equations satisfy the equilibrium equations (10), and through consideration of the 
plate constitutive equations (14), the strain-displacement relationships (5), and the pertinent 
boundary conditions, it can be shown that the associated displacements satisfy 

u L , x = = constant 

v L , y = tyi = constant (46) 

u L<y ~ u Uxx = v Lx = v Uyy = ® 

The VIPASA buckling analysis uses the unit strains e xt and £*, but does not require the 
displacements 14 and v L . The determination of the unit strains is straightforward following, for 
example, the procedure used in PASCO [2], Where plate strips join, the issue of joint equilib- 
rium boils down to the requirement that all plate strips comprising a continuous initially flat 
skin carry the same value of N yi , and all other strips carry no N y loading. Compatibility of dis- 
placements in the x-direction Is assured by the uniformity of e x within the entire panel. 
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Compatibility of the y- and z-directed displacement components at the joints is not nec- 
essarily assured by the procedure used in PASCO. The distribution of N yL among the plate 
strips within a panel is determined in PASCO with the assumptions that the distribution is 
statically determinant, and that global loads do not induce any bending of plate strips within 
the structure. In order to make these assumptions, the PASCO procedure ignores the issue 
of displacement compatibility at the joints in the y-z plane. This is acceptable for singly con- 
nected structures, but, for multiply connected structures (for example, panels with hat-shaped 
stiffeners), compatibility of the displacements at the joints will not be satisfied exactly. If dis- 
placement compatibility were properly accounted for in such structures, even a perfect struc- 
ture could exhibit slightly nonlinear behavior from the onset of load application. This 
discrepancy is noted here, but for practical structures the discrepancy may be insignificant. 
Along the same vein, if N yL loading is imposed on a panel skin which incorporates joint-line 
eccentricities (as would be the case if stiffener flanges were being modelled) or is imposed 
on the skin of an unsymmetrically stiffened panel, then bending-extension coupling would exist 
in the physical structure which would not be simulated In the prebuckling analysis Some of 
this nonlinear behavior can ultimately be captured by the nonlinear analysis. 


3.3.1 Formulae for Determining the Prebuckling Solutions 


The specific formulae for determination of the two linear, unbuckled response states are 
presented here. They are compatible with the equations used in PASCO [2J but are redevel- 
oped here for completeness in the documentation, using notation consistent with the present 
development. As stated before, biaxial loading is permitted only if there is a continuous planar 
skin connecting the boundary node lines; otherwise only loading N & is permitted. Despite the 
limitation on the permissibility of a second load component, it is convenient to develop the 
solution assuming that a planar skin exists and that biaxial loading is imposed. Once this 
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solution is obtained, various simplifications can be made to handle special cases of geometry, 
loading, and boundary conditions. 

A vector { k } of length P is used to specify which plates are part of the panel skin. The 
elements k p of vector {k} are defined such that 

{ 1 if plate strip p is part of the panel skin 

P-1.2 P (47) 

0 otherwise 

The unit global load Wxg is the mean unit axial load in the longitudinal direction per unit 
width of the panel, and can be expressed as 




(48) 


P = 1 


where B is the width of the panel between the boundary node lines, A/ft is the unit value of 
N x on plate strip p, and b p is the width of plate strip p. (In the above equation and in the fol- 
lowing equations, the superscript p is not intended as an exponent.) The unit global load Nyc 
is the unit normal load per unit length of the panel, and it acts on the panel skin at the 
boundary node lines. This load is carried by all plate strips in the panel skin, so that the unit 
y-normal stress resultant in plate strip p is given by 




(49) 


The unit normal stress resultants within each plate strip are related to the unit normal mid- 
surface strains of the plate strip through the plate constitutive equations (14): 


n x l = a u s x l + a 12 


(50a) 




(50b) 


where c xL , the unit normal strain in the longitudinal direction, is uniform throughout the panel. 
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Using equation (49) and the equation (50b), the y-normal mid-surface strain in each plate 

A 

strip can be expressed In terms of the unit longitudinal strain and the unit load Nye- 



Ny G X p -c x A 


p 

12 


A 


P 

22 


(51) 


The unit longitudinal strain associated with the specified global load components can be de- 
termined by using equation (50a) and equation (51) in equation (48) to obtain the following 
expression: 


A 



N y G S 2 

Si 


(52) 


where the constants Si and S 2 are given by 



(53) 



With the unit axial strain now known, equation (51) is used to determine the unit y-normal 
strain within each plate strip. Equations (50) are then applied to determine the unit normal 
in-plane stress resultants within each plate strip. 

On occasion, it is of interest to know the change in the width of the panel, which is the 
change of dimension between the two boundary node lines. The change in width, denoted by 
Av l , is simply the sum of the changes in width of the plate strips comprising the panel skin, 
and this can be expressed mathematically as 
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( 55 ) 


p 



p = i 


Define the mean y-normal strain of the panel skin to be 



(56) 


Using equations (51) and (55) to re-express equation (56), the mean y-normal strain in the skin 
can be expressed in terms of known values: 



(57) 


where constants S 3 and S 4 are given by 


p 

s 3 = /C k P b P ~7p 

p = 1 *22 


(58) 


S 4 = 



(59) 


At this point, all of the desired values associated with the linear unbuckled solution are 
known, but a few additional relationships are developed here for use when the boundary 
conditions are specified in different ways. By using equation (52) in equation (57), the normal 
unit load on the side boundaries can be expressed as 


A /a a S 4 \f BS ‘ | \ 

Ny ° = + N ^~s; X S^a + SaS* ) 


(60) 


From equation (52) it can be determined that 
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( 61 ) 


a N B-c S, 

N — — — — 

* S, 


Equating the above two expressions for A/*;, the following expression for A/* can be obtained: 


N x =E X 

x G x L 


Si S3 + S2S4 

BS^ 


+ £ 


So 


(62) 


For configurations having no continuous planar skin, equation (62) reduces to 


A A Si 

N = £ — — 

X L B 


(63) 


3. 3. 1.1 Boundary condition options for a panel with a continuous , planar skin 

There are four different cases identified for sets of parameters which may be used to 
specify boundary conditions for the unbuckled panel. The first three cases accommodate the 
various options for specifying boundary conditions which are to be in effect in the regime of 
nonlinear response (see Table 1 on page 42). The fourth case is for the additional prebuckling 
solution used to adjust the ratio of the global load components. The four cases are discussed 
individually below. 

Case la) N xG and N y o specified: For this case, equations (49) through (59) provide the 
solution. The sequence of application of the equations is (53), (54), (52), (51), (50a), (49), (58), 
(59), and (57). 

Case 1b) N, 0 specified and = 0: For this case, equation (60) is used to determine the 
effective unit load Nye. Next, the sequence given for Case (la) above will provide the complete 
solution. 
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Case 1c) e,L and tyc specified: First, equation (62) is used to determine the effective unit 
load NxGt then equation (60) is used to compute the effective unit load Nya . Next, the sequence 
given for Case la) above will provide the complete solution. 

Case Id) ll L = 0, Nfc specified: This solution is required only in conjunction with Case la). 
The equations which have already been presented are applied, noting that they now apply to 
strain and stress-resultant measures having the superscript A. With constants S, through S 4 
already determined, the sequence of application of equations is (57), (51), (62), (50a), and (51). 


3. 3. 1.2 Boundary-condition options for a configuration with no continuous, 
planar skin 


These configurations are assumed to be unable to withstand a nonzero value of the ap- 
plied load NyG , as such loading would induce initial bending in the plate strips which is not 
accounted for in the current method. Indicator k p of equation (47) is zero for all plate strips. 
Only one case is identified for specification of the prebuckling load state. 

Case 2a) N xe specified, Nyo = 0: The sequence of application of equations for this case is 
basically the same as for Case (la) for a panel with a skin, except that there are some omis- 
sions. Constants Sz, S 3 , and S 4 of equations (54), (58), and (59), respectively, are all zero. The 
sequence of application of the equations is (53), (52), (51), (50a), and (49). 


3.4 Buckling Eigensolutions: 


The equations which determine the buckling eigenfunctions are obtained by setting to 
zero the terms of the expanded equilibrium equations (44) which are of first order in the modal 
amplitudes. A set of equation trios results, one trio corresponding to each coefficient q h The 
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i th eigenfunction satisfies the trio of equations corresponding to q t , with h (the I th buckling 
eigenvalue) substituted for 


NxfX + ^xyyy ■*" ^i^y L u hyy ~ ® 

^xy,<x "h Ny t ,y + h flx^bxx = ® 

Mxyxx "h 2M X y.,xy + My.,yy + X t (N^W vxx + NyWj,yy) — 0 


(64a) 

(64b) 

(64c) 


In general A, will be negative, corresponding to a compressive end loading. An infinite number 
of eigensolutions exists for the eigenvalue problem, and it is assumed here that 
eigensolutions are ordered according to increasing negative magnitude of the eigenvalues: 


^ ^2 A3 > ... 


(65) 


It is noted that equations (64a) differs from the corresponding equation used in the 
VIPASA analysis [5], the latter equation being 

NxyX ~i~ ^xyyy “h ^i^x L U i’XX = ® 

The term W xt u„„ in the above equation was obtained through the use of a nonlinear theory of 
elasticity in which N x is considered to be the primary loading in each plate strip. It Is argued 
here that the third left-hand-side term of both equation (66) and equation (64a) can be neg- 
lected in the current application. The term W xt u ( ,„ in the VIPASA equation is felt to be of 
negligible importance, based on the consideration of the effects of moderately large in-plane 
rotations which gave rise to the plate theory used here (see Appendix A). Buckling analyses 
of stiffened panels were performed using VIPASA both with and without the last left-hand-side 
term of equation (66) (involving the alteration of a single line of VIPASA source code), and the 
resulting differences in the computed buckling eigenvalues were on the order of one tenth of 
one percent. The term NyM.yy of equation (64a) is considered to be of negligible importance, 
because for a plate strip in the skin of a panel, in-plane rotations will remain small, making 
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the nonlinear plate theory of von Karman valid, thus eliminating the term, whereas for a plate 
strip not lying in the skin of a panel, N, L will be zero. The VIPASA buckling solutions can thus 
be assumed to closely satisfy the buckling equations {64) for the configurations under con- 
sideration, whether or not the term N„ L ^, XX is present in the VIPASA analysis. In the present 
application, the term will be omitted from equation ( 66 ) during applications of VIPASA, 

and the stress resultant contributions N, L and A/#, will be omitted in evaluating the term 
(N y u,y),y in the first of the equilibrium equations (10). 

The procedure used in VIPASA for determining the eigensolutions is described in detail 
by Wittrick and Williams [5], Only the details of the VIPASA procedure which are pertinent to 
the present method are discussed here. It is noted that VIPASA is capable of determining with 
certainty the fundamental buckling eigenvalue for a given longitudinal halfwave number, and 
any desired number of additional eigenvalues in sequential order (increasingly negative 
value, for the conventions used here). 

The i* eigenfunction has the following assumed form on each plate strip (where a phase 
shift in the x-direction has been applied relative to the conventions of VIPASA in order to lo- 
cate the x-domain of the panel in the interval [0, L]): 

fU/-j f £f(y) cos fmrx/L'| 

{u,} = < v , 1 ^ sin mnxjL V (67) 

liv ( J l<£,(y) sin mnxjL ) 

where m is the number of halfwaves along the length of the panel for the /* eigensolution. 
When the above displacement forms are inserted into the buckling equations, there result two 
coupled homogeneous linear second order ordinary differential equations in the functions 
Zi(y) and »j,(y), and one homogeneous linear fourth order ordinary differential equation in the 
function 4>,(y). These ordinary differential equations are presented in Appendix B in equations 
(B21) and (B9). 

The solutions {(*;,} = [£, > 7 , <£,] T to the ordinary differential equations presented in Appendix 
B have the following general form: 
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( 69 ) 


where there is, in addition, a compatibility condition between & and rj,, expressed as 


{£,} = lR jk ]{G k ) 


where [fy] is a matrix of constants, and C h E h and G, are arbitrary constants. The values of 
the constants R jk and the forms of the functions fj{y) and g y (y) all depend on h and m , and are 
tabulated in Appendix B. There are eight independent arbitrary constants, which take on 
specific values with the specification of the eight generalized side-edge displacement ampli- 
tudes, {£,e} (e = 1, 2), defined by 



UYe) . 

U(y e ) 

U(y e ) 

4>,'(y e ) 


(70) 


The generalized displacement amplitudes {& e } at the side edge of a plate strip, when 
referred to the associated node-line number n and the global coordinate axes (see Section 
3. 1.2.2), are denoted by {Or} = [Or Vp Wp ^ 7 ] . In returning a computed eigensolution, VIPASA 
returns the set of amplitudes {0,}" corresponding to each node line. (Actually, as a result of 
the use of complex variables, VIPASA returns the vector [ —Up V? Wp ^ 7 ] .) With the /-de- 
pendence of the displacement shape functions known, the eccentricity transformation of 
equations (20) can be evaluated, and a transformation can be specified relating {<*,>} and 
{Op} [5]: 


ile) = lT ecc ]lT r l{U?} 
6 


(71) 


where n is the node line number corresponding to edge e of the plate strip under consider- 
ation. [T r ] is the rotation-transformation matrix given by 
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( 72 ) 


[^ = 


1 

0 

0 

0 


0 0 0 
cos p — sin 0 
sin \x cos fi 0 
0 0 1 


and [Tecc/] e 1$ the eccentricity-transformation matrix for the buckling eigensolution, given by 



1 

0 

0 

0 


mn „ 
f_ e ye 


1 

0 

0 



(73) 


In order to assure that the buckling eigensolutions satisfy side-edge boundary conditions 
and equilibrium of the node-lines (the plate free edges and side-edge joint lines), a transfor- 
mation of the type just presented for the generalized edge displacements must be available 
for transforming the generalized force-resultants at the edge. The generalized force-resultants 
of equations (13), {4} = [f xe f ye ht m e ] r , are expanded in the same manner as was done for the 
equilibrium equations, to get the portion of {f e } associated with the I th buckling eigensolution. 
The portion is denoted {/,*}, and is given by 


{fie) - 


n Ay, I 

n yNy t \ 


y e 


y. 


n y (2M X y, X + My r y + XjNy^ by ) I . 


Ye 


- HyMy. I 


These can also be expressed in terms of separated variables: 




; ___ mirx 

L COS — ; 

A i9 

r} ein rrrnx 

I f y,e s,n ~[~ 

I ; mnx 

K sm ~r 

mnx 


m ie sin- L 


(74) 


(75) 
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where the amplitudes {}«} are scalars which depend on the functions (6(y)}- When the ampli- 
tudes { f, e } are referred to the node line and global coordinate axes, they are denoted {f, p }, and 
the x-dependence associated with the latter quantities is the same as for the former quantities 
(equations (75)). For {fr} to be statically equivalent to { f ie }, the virtual work of {f ie } acting 
through virtual displacements { 6u, e } will be equal to the virtual work of {fr} acting through 
virtual displacements {e5Lff}. This condition is met by requiring that [5] 

vti=ir,i r n.cc?y») ™ 


The associated generalized node-line force-resultant amplitudes, {Fr}, are equal to the sum 
of the contributions {fr} from all plate strips terminating at a node-line, as expressed in 
equation (24). The values {Ff} are zero at all nonboundary node-lines (equation (21)). At 
boundary node-lines, {Ur} and {Ff} satisfy the homogeneous forms of the boundary conditions 
selected from Table 1 on page 42. 


3.5 Second-Order Displacement Fields 


The equations used to determine the functions {uy} = [u, ; v,y w v ] r for each plate strip are 
obtained by setting to zero the terms of the expanded equilibrium equations (44) which are 
of second order in the modal amplitudes, excluding the terms which contain the parameter 
Ay. (The excluded terms have a different characteristic waveform in the longitudinal direction 
than is used for the shape function {u, y }.) Each resulting ij* trio of partial differential equations 
is to be satisfied independently, and these are given by 

^Xjj'X + y.^y + “ 2 ~ [(flyUj t y),y + (NyUj,y) iy ] — 0 {77 a ) 
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^xyij'x + Ny.j.y + ^t>^x L v i]'xx "l" 2 [( N xYj<x)<x (NxVi<x)<x^ ® 


(77 b) 


M Xj j'XX + 2/W xv , y + My.yyy + X b ( W^VVy,^ + NyWy.yy) 

+ ~2 t-(^Xi W j>x ^xy t w j'y)’X (^xy^j'X + ^y^ py)-y 

+ (Nxj W i<x + ^xy^fy)’X + (^xy^i<x + ^y / w 'i<y)'y3 = 0 


(77c) 


where the second-order contributions in equations (44) were regrouped so as to obtain 
equations which are symmetric with respect to the subscripts / andy. The term XbN y iUtj iy y has 
been omitted in the first equation above, consistent with the discussion in the Section 3.4 . 
The parameter X b is a reference value of X, not necessarily the exact value of X at which an 
equilibrium solution is sought. The selection of the value for X b will be discussed in Section 
3.7 . The above equations involve known eigenfunctions {14} and { 14 }, and the unknown 
shape-function {u, y }. Let m and n denote the number of halfwaves in the x-direction for 
eigenfunctions {u,} and {u y }, respectively. Separation of variables can be used to convert the 
trio of partial differential equations (77) into two trios of ordinary differential equations in the 
variable y, by assuming the following form for {t 4 y } [56]: 



i a iy) sin mnxIL 
*1 (y) cos mnx\L 
<f> a/ (y) cos mnxjL 


(78) 


where 


{m + n a = 1 
m — n a = 2 


(79) 


The ordinary differential equations which result are derived in the following sections, along 
with various other expressions which must be known for a full determination of the second- 
order displacement fields. First the generalized strain measures and the stress resultants will 
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be expressed in terms of separated variables, and then the governing equations will be ex- 
pressed in terms of the unknown displacement shape functions, again with separated vari- 
ables. The generalized force-resultants at the side edges will be expressed in the same 
fashion, as needed for the imposition of boundary conditions. 


3.5.1 First-Order Fields in terms of Separated Variables 


The first-order displacement fields (the buckling eigenfunctions) have the following form, 
copied from equation (67): 

r £/(y) cos rmrx/L'j 

{u,} = < v, > = < rti(y) sin mnX l L 7 ( 80 ) 

(.w,J l<£,(y) sin mnxjLJ 

where the functions {<*,} = <f>,Y are known based on results from a buckling analysis. The 

first-order mid-surface strain and curvature fields are given in equations (37) and (38), and are 
copied here: 


1 

r ii. > 

| | 

1 | 

r w i«x ) 


- u i’y + v i<x-‘ 

[ W H 

1 w l 

1 w fyy ( 

L2w„ xy J 


(81) 


Using equations (81), the mid-surface strain and curvature fields can be evaluated and ex- 
pressed using separation of variables: 


{ e x .(y) sin mnx/L x 
tyfy) sin mnx/L 
Y xy ly) cos mnxIU 


{ ic x .(y ) sin mux/L -j 
Ky.iy) sin mnxjL > 
1c xy.(y) cos mnxjLJ 


(82) 


where 
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E *, = - m j; it 
i v =n( 

yxy, = £/ + m *u 


(83) 


K *! = m2 (f) 4 1 

*y = -4>i" 

*xy, " “ 2 m •f" <£/' 

Now the first-order stress resultant fields can be expressed using separation of variables: 


N x .(y) sin mnx/L 

,M x fy) sin mux/L 


{N,} = -j N y [y) sin mnxjL J> 

{M,} = < M y (y) sin mnx/L l 

(84) 

^N xyi (y) cos mnxlL J 

^M xy .(y) cos mnx/L J 



where the y-dependent contributions to the stress resultants depend on the y-dependent 
portions of the strain fields: 

{■ N ,} = [4]{e,} {M,} = [D]{k,} (85) 

where the plate stiffness matrices [4] and [D] are given in equations (15). 


3.5.2 Second-Order Fields in terms of Separated Variables 


By substituting the separated form of the unknown functions {uy} from equation (78) and 
the known eigenfunctions { u ,} from equation (80) into the second-order mid-surface strain and 
curvature shape-functions (repeated here from equations (37) and (38)), 


r U/y.x+1/2 {Vi, x v j , x +w i , x w j , x ) -j 

{«,;} = < v ij~y + 1 /2(U/,yU,.y + W,, y Wj. y ) V 
t U/y.y + Vij'X + 1/2(w f , x w y , y + W hy Wj, x ) J 


{*//} = - 


(86) 


and applying the following trigonometric Identities, 
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(87) 


2 

sinmC sin nC = 1/2^ - s a cos mC 
** 1 
2 

cos mC cos nC = 1/2^ cos mC 

a= 1 


sinmC cos nC = 1/2^7 sinmC 
0= 1 
2 

cos mC sin nC = 1/2 ^s a sinmC 


oc= 1 


where s. is given by 


f 1 if a = 1 

j-1 if a = 2 


( 88 ) 


the strain and curvature functions can be expressed in terms of separated variables as fol- 
lows: 


2 

COS '”* X / L 1 

2 rK xo ..(y) cos mTix/L -j 


(«#>-£- 

S ya .(y) cos mrrX/L > 

{*,;} = X j ^ya/y^) cos ™ rtX / L [ 

(89) 

«= 1 

'xy«tj(y) sin m^rx/L-' 

“ =1 ^xya tJ (y) Sin mnxjU 



where the y-dependent portions are given by 


= + mn { "L ) + W 

^=< + T ( W" s «W) 

W,, = V + ( m W + 


^ = ^ 2 (f ) K n 

K yctjj = ~ ^tXjj 


(90) 


Now the second-order stress resultant fields can be expressed using separation of vari- 
ables: 


0 

N xa .( y ) cos mnx/L 

/Vf x (y) cos mnx/L 
2 r 11 

II 

M r 

W yaj/ (y) cos mnx/L l 

{M tJ } = /W ya ..(y) cos mnx/L 

or= 1 

N xya .(y) sin mnx/L 

M xya (y) sin mnx/L 


(91) 
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where the y-dependent contributions to the stress resultants depend on the y-dependent 
portions of the strain fields: 


{Ay = [*]&.} {M aij } = lDl{K aij } 


(92) 


where the matrices [4] and [D] are given in equations (15). 


3.5.3 Governing Equations in terms of Separated Variables 


The displacement and stress resultant shape-functions appearing in the partial differen- 
tial equations (77) are now replaced with the corresponding expressions given in equations 
(78), (80), (84), and (91). Through use of the trigonometric identities of equations (87), the par* 
tia! differential equations can be expressed in terms of separated variables: 


2 

2 { - ™ f + W xya/ ; + 1 IS'lNyfiy + (Ny, */)']} sin m*x/L = 0 
«= 1 


(93) 


5 "! \ m l_ N xy<*ij + ^y«ij *b^x L m ( 

a= 1 U 

+ -j ( ) [m(n + s a m)N x ,> J( + n(m + s a n)N x ^]j cos mrrx/L = 0 


(94) 


Xi | - 17)2 { ~L ) M xa ii + 2m T" + + _ Wjc r m ( ~L ) ^ a ‘i + N ^ a 'i 

+ j^m(s a m + n)( -J- ) N Xy <£, + (2m + s a n) -J- N xy , </>,' + m~ N xyj '4>i - s a (W y . <£,')' 

+ n(m + s a n)( j- ) N x <f> i + ( s a m + 2 n) j- N xy .4>/ + n y- N xy ’4>, - s a (N y .4>j'Y j 


x cos mnx/L 


= 0 


(95) 
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In each of the three equations above, the two x-dependent trigonometric functions that 
correspond to « = 1 and a = 2, respectively, vary independently and are in general nonzero. 
Thus the two y-dependent portions of each equation, corresponding to a = 1 and a = 2, re- 
spectively, are independently set equal to zero. This provides six ordinary differential 
equations in the variable y for each index pair (//). (The one exception occurs when m = n 
making m - 0 when a = 2, so that sin mnxIL = 0 for a = 2 in equation (93), In this case there 
are only five equations obtained.) 

To get the final working form of the six equations, the stress resultant terms in equations 
(93)-(95) are expressed in terms of {{,} and {£„ v } using the equations developed in this section. 
This provides two uncoupled trios of equations governing the functions {<*.,.,} (a = 1, 2). The 
equations are linear, nonhomogeneous ordinary differential equations with constant coeffi- 
cients, having the following form, where the ij subscripts are omitted: 


Cut«" + C t 't a + C u r, a ' = F a (y) 


(96a) 

DuS a ' + W + D uVtl = G a (y) 

« = 1, 2 

(96b) 



(96c) 


where the constant coefficients are given by 
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W = F a ({ti(yMZj<y)}) 

= ± \mnm( f f A u (^j + + mf *««,'«/ - sjtf’) 

L ( 98 ) 

- s a j- A 6e l{m + s a ri)4> {<$>)' + mW + s a n<£ "<£,] + (m + s a n) -J- -4 12 £,'{/ 

+ f 4 12 (ro<tf/' + s a nZ,"Zj) - 4 22 [( n/ "£/ + n/«/0 + «.«/'•»/ + «/•»/')]} 


G a (y) = G»({iM.{p) 

{ 2 2 

“ V”( Y ) *66 ( m <t>i<t>/ + s a n<t>i<t>j) ~ mn( y- ) /A 12 (^/^ + »?,»)/ + 4>!4> t + 4>,4>j) 

- A 22 i(trz; + tit") - Sait,"*/ + wpi ( 

+ mn( -J- j A^(m£it]j + rnj^y) - mn| — J + ij,»j/) 

+ S a mn ( ^ ) ^ 1 l( n ^ + m »J^)-S a (-^) ^ 12 ( n Vny + m\»/y')|- 


w a (y) = w a ({^(y)}.{^(y)}) 

= - y jn(m + s a n)( j- ) ( - m j- A^Zrfj + A^^j) 

+ (s a m + 2 n) f A^ttW + m f Vitj) + nf A 6 tf,"*] + m f 'W 

+ s a m f A n {Z’4>; + {,*/') - s a A 22 hr<t>j + Vi'4>j") (100) 

+ m(s a m + n)( j- ) ( - n -j- A^fi + 4 12 <M/) 

+ (2m + s a n) f Adfifti + nf + mf A m (<j>{/’ + nf <M/) 

+ s a n f A u Wlf + Wtj) ~ ,/)} 

In equations (97)-(100) both m and s. depend on the value of a. For the special case where 
m = n and a = 2, equation (96a) is omitted and (^s0. 


Method of Analysis 


7! 



3.5.4 Generalized Displacements and Generalized Force-Resultants 


Along Boundaries: 


Expressions must be developed for the generalized edge force-resultants associated with 
the second-order displacement fields, to facilitate specification of the boundary conditions 
which accompany the differential equations (96). The portions of the generalized force- 
resultants at the side edges of a plate strip ({/»} of equation (13)) associated with the second- 
order displacement fields are given by 


~ ( 1) ^Xy,y + 2 (NyU),y+ Ny.Uj,y)] 

f Zjjl = ( — 1)®[2 M xy .,, x + + ~ *f NyWj,y + N xy Wj, x 

^y e = -(-1)X,l 

e = 1, 2 


+ NyWj.y ) ] 


y. 


( 101 ) 


The above equations can be expressed in terms of separated variables by substituting in for 
the stress resultant and displacement terms the separated-variable expressions presented 
earlier in this section, and applying the trigonometric identities given in equations (87). 
Equations (101) can thus be converted to the following form: 


L ljt sin ™ nx i L 



m. . cos mnxjL 


where the scalar amplitudes introduced above are given by 


(102) 
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L, f = ( -V e lN xya ^ + -J- (s a N yj £ ’ + N y t/)1 1 

■^e 

4a ,y< = ( — ^ xy*ij "*■ My a/ j + ^bNy^ajj 

+ ^ C ( m N X y4> I + nN X y,(f>j) — S a (Ny.<f>l + N y<f>j )]} 


( 103 ) 


To get the final working form for the above equations, the stress resultant terms are expressed 
in terms of {£,} and {£ >(7 } using the equations developed in this section. Omitting the ij sub- 
scripts, the equations have the form 


L = (-i) e Cc 1 ^ a ' + c 2t n a +F a (y)3| 

L = ( -V*l\ t a + D 2s + G a (y)] | 

y ' a = 1,2 (104) 

L = ( 4>a'" + *2. *«' + «.(/)] | 

m a , = -(-1) e [E 3 > a '' + M ] l 

where the constant coefficients are given by 


Ci« — ^66 

— A 

C u = ~ m 

f A 66 

— A n 

D u — m — A 12 

®2. ~ A 22 


E u = — D 2 2 

E,. - m 2 ( 

2 

f) (d 12 + 4D K ) + ;yj yi 

= - d 22 


f)‘»- 


and the functions of y are given by 


(105) 
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F a (y) = M.tf A») 

= -J- [s a f Admit*/ + sMi'tj) ~ f AnWitf + + *22(n/«/ + *« W)] 

G>) = G a ({^(y)},{^(y)}) 

= _ 4'[ mn (f') Awt’li’lj + ‘f’i'f’j) + A 22(Zi't/ ~ Sa^/*^/) j (I®*) 

H a (y) = « a ({^(y)}.{^(y)}) 

= - J- [ f - A 66 (n(,'<i>j + m </>,{/) + mn ( f - ) AseMy + <My) 

+ s a -£- A 12 (m^i<f>j + n<j>/£j) — s a A 22 (tji'4>j' + </>/'»?/)] 


For the second-order displacement fields, the four components of the associated gener- 
alized displacements and generalized force-resultants along the plate side-edges and the 
node lines follow a consistent pattern with regard to x-dependence, shown in the following 
equation, where ij subscripts are omitted: 


{u e U n f e f n F 


(L U” f xa . & F Xa ) sin mnxjL 
„ _ y[(l V* fy» i Fy a ) cos mnxjL j 
} «= i|(<L f zae f z n a F Za ) cos mnxjL J 

(k, M*) cos ™ nX l L 


(107) 


where the superscript n denotes the number for the node line corresponding to side-edge 
number e of the plate strip under consideration, and ( — ) denotes scalar values. By applying 
the transformation equations (19) and (20) in the same manner as was done for the first-order 
displacement fields (the buckling eigenfunctions), the rotation and eccentricity transformations 
can be expressed as matrix operations involving the scalar amplitudes appearing in the 
right-hand side of equation (107). For specific values of ij, and a, the transformation of gen- 
eralized displacements has the form 

= LT, cc ]JLT r -]{U n } e = 1,2 (108) 
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where the subscript a is now also omitted, the rotation-transformation matrix [T f ] is given in 
equation (72), and the eccentricity transformation matrix is given by 


1 

0 



0 


A TC A IT « 

m j_ e ye m j_ e ze 0 
1 0 e ze 

0 1 - e ye 

0 0 1 


(109) 


where e ye and e ie are the eccentricity measures (defined in Figure 4) for edge e of the plate- 
strip edge under consideration. Because of the presence of m in equation (109), the 
eccentricity-transformation matrix depends on the value of a. The matrix is different in form 
from the corresponding matrix which applies to the buckling eigenfunctions (equation (73)). 

Using the same approach as was used in Section 3.4 for the buckling eigensolutions, the 
generalized edge force-resultants are transformed from the side-edge of a plate to the asso- 
ciated node line and global coordinate axes by using the following equation: 


{h = wT wC&) 


( 110 ) 


The amplitude of the generalized node-line force-resultants, (P), is the sum of the contrib- 
utions { } due to the plate strips terminating at the node line (see equation (21) ). 


3.5.5 Boundary Conditions Along the Node Lines 


Along nonboundary node lines it is required, consistent with equation (24), that {P,*J} = 0. 
This is enforced by imposing the condition {p, v } = 0 (a = 1, 2), where {F;,j} was defined in the 
preceding paragraph. Along boundary node-lines, it is required to satisfy the homogeneous 
form of the boundary conditions which have been selected from Table 1 on page 42, with the 
exception of the third boundary condition option for displacement component V n , For this op- 
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tion it is required that V (> ,! = 0 and F^- 0. It can be shown that these conditions will be sat- 
isfied by requiring 


V" = 0 if a = 1 

v" =0 ifm^rO ( 111 ) 

L = 2 

Fy.,j = 0 if m = nJ 

3.5.6 Solution of the Boundary Value Problem for {£„(/)} 

Because of the complexity of the expressions found in \.F tll iy) G tlJ (y) H a// (y)] of equations 
(98)-(100), it was opted to perform a finite difference analysis to obtain values of the functions 
{£«,;(y)} at a finite number of evenly spaced y-staiions on each plate strip. The finite-difference 
analysis procedure is presented in detail in Appendix C. 

3.5.7 Symmetry of Second-Order Displacement Fields with respect to i 
and j 


The boundary-value problem for the displacement fields {u, v }, /,/* 1, 2, .... has been 
posed so that {u /y } is symmetric in / and or expressed symbolically, 

{«#}-{«//> < 112 ) 

This raises the question of how the functions {£ a ,y} and {£*/,} relate. The expression for { 11 $ 
of equation (78) can be used to re-express equation (112): 
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(113) 


LM sin(m + s a n)rrx/L 


( a ..(y) sin(n + s a m)nx/L 


2 ' ) 2 f p “ 

X ] cos ( m + V’) 7rX / L f = Z ) cos(n + S « m ) 7rx / L 

a= 1 C J a= 1 L 


4> a/y (y) cos(m + s a n)*x/L 


4> a ..(y) cos (n + s a m)nxjL 


where the following substitution has been made, based on equations (79) and (88): 


m = (m + s a n ) 


(114) 


It is noted that 


(m + s a n) = s a (n + s a m) 


(115) 


and that 


sin[s a (n + s a m)ffx/L] = s a sin[(n + s a m)?rx/L] 
cos [s a (n + s a m)n-x/L] = cos[(n + s a m)?rx/L] 


(116) 


Using equations (115) and (116) in conjunction with equation (113), it is easily shown that 


^,(y) = ^(y) « = i 

« = 2 

%(/) = %(/) ] 

" " V a = 1,2 


(117) 


Computational efficiency is gained by avoiding redundant computation of {uy,} if {u, ,} is known. 
For example, to get the full set of required functions, compute {try,} only if j > i. 
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3.5.8 Discussion of the Load-Dependence of the Second-Order Fields: 


Three special cases are identified which effect the load-dependence of the second-order 
displacement fields. The first case is for an unstiffened panel configuration which buckles 
out-of-plane. In this case, the only nonzero eigenfunction component is w>. This causes H.{y) 
of equations (96c) and (100) to be identically zero, making the differential equation (96c) ho- 
mogeneous and accepting of the solution <j>^(y) = 0. Furthermore, the conventional von 
Karman nonlinear plate equations can be assumed to apply, with the consequence that the 
coefficient D u of equation (97) lacks the term containing the load parameter, and the equations 
governing and >j a , v become load-independent. 

The second case is a stiffened panel undergoing local postbuckling, such that in the 
buckling eigenfunction only the component w, will be nonzero on each plate strip. In-plane 
rotations of each plate strip will be small, making the von Karman nonlinear plate equations 
acceptable. As was the case for unstiffened panels, discussed above, the differential 
equations governing and tj, t j (equations (96a) and (96b)) are load-independent, and the 
differential equation governing (equation (96c)) is now homogeneous. However, unlike the 
case of the unstiffened panel, the three equations (96) are coupled through the displacement 
compatibility and equilibrium conditions at the joints between noncoplanar plate strips. Thus, 
the second-order fields have three nonzero components, and because the reference load pa- 
rameter At, appears with equation (96c), the second-order fields are load-dependent. However, 
if the "classical assumptions" for local postbuckling analysis (see Section 2. 2. 2. 3) are invoked, 
then the coupling of equations (96a) and (96b) with equation (96c) is removed, and the 
second-order displacement fields become load-independent, with component <£«,, equal to 
zero. The "classical assumptions" referred to are not incorporated in the analytical results 
presented in this document, but as will be discussed in Chapter 4, when performing local 
postbuckling analysis, the assumption is made that the second-order displacement fields are 
load-independent. 
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The third case is when m = 0, which occurs when m — n and « = 2 (see equation (79)). In 
this situation, as stated earlier, £./y = 0. Furthermore, equations (96b) and (96c) reduce to 


^22^a" ~ G a (y) 

- D 22 4> aii "" + hNy L = H x (y) 


(118) 


and thus, if N yL is zero (uniaxial rather than biaxial loading) then the displacement contrib- 
utions tj./j and 4>iij are load-independent. 


3.6 Enforcement of the Load Ratio 


This discussion applies to cases in which the structural configuration admits biaxial 
loading and the unit global load Nyo is specified. For this situation the ratio of the global load 
components is to be maintained at a constant ratio, as specified in equations (23) and re- 
peated here: 

AAA , . 

(119) 

In the linear regime of response of the perfect structure, the equilibrium solution which also 
satisfies equation (119) is obtained by setting the load parameters and the modal amplitudes 
g, (/ = i t 2, ...) to zero in the expression for displacements of equation (33). In the regime of 
nonlinear response, however, the load parameter A is needed to assure that equation (119) 
is satisfied. 

Assume for the simplicity of this discussion that the panel skin terminates as a simple 
plate edge at the two boundaries node-lines, so that the normal force-resultant at a side-edge 
of the panel is simply f y = nyN y , evaluated at the appropriate location. Global force resultant 
Nyc is evaluated at a boundary node line, and it has the definition 
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( 120 ) 



where n i is the node-line index number of a boundary node-line. Global force resultant A/«g is 
evaluated at x = 0, and it has the definition 


N ^ = 



( 121 ) 


where b is the width of the p* plate strip, B is the reference width for the panel, p is the 
plate-strip index number (not intended as an exponent), and P is the total number of plate 
strips. 

By expressing N x and N y , appearing in equation (121) and equation (120), respectively, 
with the expansions given in equation (39), the global loads A/*; and NyG can be expressed in 
expanded form: 



where N x G , Nyc, and W& are all specified unit loads, A/jfc is known from an earlier analysis 
(Section 3.3), and 



(123) 
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- -1 f L % I 
L Jo AL 


It is found upon evaluation of equations . (123) that because of the form of the harmonic x-de- 
pendence of the stress resultants, 


N^hO /- 1,2,3,... 


N x = 0 U, *-1,2,3. ... 

G ijk 


The expanded expressions for the global force-resultants of equation (122) are now used to 
express constancy of the load ratio (equation (119)) to get 


HNyC - ™ xC ) + (A- A °)(Ny 0 - RN? g ) + (q, - qf)N yQj + (q flj - qfqfWy^ ~ 
+ (q,q l q k -q°q°q°)N yG ^ + - qfqf q k q^)(N yQ ^ - RN^J = 0 


The above equation will be satisfied if A is taken to be 


A = q,A + qfyAy + q,qjq k A l]k + qflfl^A tjk{ ij, k, i = 1, 2, 3, 


where 


A A A A 

AP — RN 


n: - RN : 


Ny G - RN- 

yQ ij // 

A J» A x 

AP — PAP 

' v y G 


w y G 

yG >jk< ^ijkt 

A A a X 

AP — PA/ 

*g 
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Now the variable parameter A can be replaced wherever it appears, using the form given in 
equation (128). 

One unanticipated use for the parameter A was discovered through applications to test 
cases. It was originally believed that the need for the parameter was limited to cases for which 
a nonzero secondary load component Nyc Is present [62], This belief was based on the fact 
that the shape functions { u L } satisfy the boundary conditions, and the shape functions { u ,} and 
{u,j} each satisfy the homogeneous forms of the boundary conditions for their respective 
boundary value problems. However, the boundary-value problem governing the 
eigensolutions { u ,} actually governs the buckling of an infinitely long structure having periodic 
supports, rather than finite-length panels which may be of interest. Consider a wide (relative 
to the supported length) flat panel, stiffened on one side with blade stiffeners. During global 
buckling of an infinite-length, periodically supported panel, the skin will have an incremental 
longitudinal strain which alternates between positive and negative values in sequential bays 
along the loading axis, due to the eccentricity of the skin relative to the neutral bending axis. 
Away from the side-edges of the panel, this alternating longitudinal strain induces a normal 
stress resultant N y which alternates between positive and negative values In sequential bays 
along the loading axis. Such a panel would typically be represented in an analysis by a single 
stiffener attached to skin which extends to either side the distance of one half of the stiffener 
spacing, with symmetry conditions imposed at the two side-edges. The global buckling mode 
for this model exhibits nonzero normal stress resultants at each side boundary, which provide 
zero net side loading on each boundary when integrated over an even number of bay lengths 
in the longitudinal direction, but which provides a nonzero side load when integrated over only 
a single bay length. Using the current method to analyze this representative width of a 
finite-length panel, the parameter A of equation (129) is nonzero, acting to modify the buckling 
mode-shape to ensure that the average load N y over the length of the the single panel bay is 
zero. 
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3.7 Final Determination of Equilibrium Solutions 


With the shape functions determined as described in the previous sections, it remains to 
determine final, approximate equilibrium solutions. The originally planned method for ac- 
complishing this was unsatisfactory, and it is of some interest to note the sources of detri- 
mental error in the original approach. For this reason, the original approach will be discussed 
briefly, followed by a discussion of the final approach used in obtaining results. 


3.7.1 Original Approach 


In the original approach to obtaining final equilibrium solutions, M eigensolutions are first 
selected to establish finite basis for the expression of displacements (see equation (33) ), 
strains, stress resultants, and so forth. The next step is to form a set of M nonlinear algebraic 
equations through use of the following weighted-residual expressions: 

p 

/* ' [{EWk'X + W X yy "F (WyM.y)ty]U/ + I- N x y, x + Ny,y + (W x v, x ) ljt ]v ( - 

P - V lM X , XX + 2 M X y, X y + My,yy + (N„W, X + N X yW,y), X + (N Xy W, X + NyW,y),y] W j) d A\ ^ (1 30) 

= 0 , i = 1, 2, ... M 

where the bracketed expressions now represent the residuals associated with the finite-basis 
expression of the three Euler equations (10). The weight functions are seen to be the three 
components of the buckling eigenfunctions. For a simple plate problem, equation (130) de- 
generates to an application of the Galerkin weighted-residuals method, and is suitable for 
determining exact series solutions. This is not the case for linked-plate configurations. By 
comparing the above expression to the second expression of the virtual work principle, given 
in equation (28), the source of errors can be illuminated. 
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First, a proper expression for the components of [Su] appearing in equation (28) would 


be 


M / M \ 

{6u} = Y m M + Y 2q t u J ) (131) 

/= i V j = i / 

whereas the use of equation (130) neglects the contribution of the second term in the paren- 
thesis above. A second approach was investigated in which the three eigenfunction compo- 
nents used as weight functions in equation (30) were replaced with the corresponding 
components of the expressions inside the parentheses in equations (131). However this ap- 
proach was still unsatisfactory for some problems. 

The second source of error comes in neglecting the boundary integrals present in 
equation (28). In solving the boundary value problems governing the linear, buckling, and 
second-order displacement solutions, respectively, the associated force-related boundary 
conditions were satisfied exactly. However, in computing higher-order contributions to force 
resultants at the panel boundaries and internal node-lines, there are finite errors in the sat- 
isfaction of force-related boundary conditions of equations (30) and (31), such that the errors 
must be accounted for in expressing equation (28) in order to avoid gross errors in the final 
solutions. The above-mentioned sources of error were both discovered and confirmed through 
application to a variety of test cases. 


3.7.2 Final Approach 


The final approach used to obtain approximate equilibrium solutions involves a direct 
application of the virtual work statement given in equation (32). A set of M eigensolutions is 
selected to establish a finite basis for the expression of displacements, generalized strains, 
stress resultants, and so forth. Some strategies for selecting the set of M eigenfunctions are 
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discussed in the Chapter 4. The virtual work statement of equation (32) is evaluated with the 
expanded forms for the generalized strains and stress resultants which have been developed, 
and expressed in the following form: 

M 

6W=Y J K i Sc ’i = 0 ( 132 ) 

/= 1 

where K, is the expression forming the coefficient for <5q,. Because of the independence of the 
M virtual modal amplitudes <5g„ M independent equations are obtained from equation (132): 

K, = 0 , / = 1, 2, .... M (133) 

The set of M equations (133) can be expressed as a set of nonlinear algebraic equations in the 
following form: 


M 


M M 


M M M 


xc t~ Y q j° c i ~ Y Y Y Y 

}=y J = UbI ]m 1* = V=1 


M M 


7 = 1 
M M 

7=1K = 1 
M M M 


+ X ^1 Cf J + _ X ~ X! Yj q ° q s c v 

k = 1 * = 1<f = 

M 

C iJk + ~ ^ fy^ljk 

+ ^ ^ + °to/) = o , i = 1, 2, ... 


M 


7 = 1* = i/«i 


(134) 


where the sub- and super-scripted coefficients C are constants. Before giving the expressions 
which define the constant coefficients, some shorthand notation is introduced. First, it is noted 
that, because {u, ; } = {u ;V }, the definition of the strain shape-functions given in equations (37) 
assures the following subscript reversibility: 
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{ £ i;} = { £ yi) 

( E ijk) = { E ikji 

{ E ijktf} = i c ikfk) = i E kflj} ~ l c kfjl) 


( 135 ) 


The same subscript reversibility applies to {k,,} of equation (38), {A/,y}, {A/,,*}, and { N ijkt } of 
equation (39), {M,,) of equation (40), and to Ay. Ay*, and Ay*/ of equation (128). A tilde over 
certain variables is used to denote the following abbreviations: 


IN,,} = {N,j} + A.jiN*} (136) 

{N,jk} = lN iik } + A llk {N*} 

{7,} = {c,} + A{ £ ?} 

&j } “ ty) + A >M^ (13?) 

3 ( £ p) = ({ £ /y*} + Ay*( £ L)) + 2({ E *,;} + Awy( E ?}) 

Afiijkt) = 2({ £ /y7(^} + Ay*/{ £ (.}) + 2({ E ktij} + A k f,i{c L }) 


V” = V" + A,vf 

v?i = v? j +A IJ v£ n 


(138) 


3 Ay* = Ay* + 2 ^*/y 

AAjjM = 2 A y w> + 1A ktij 


N y s r% + A ^y a 

% = % + Vy c 

%r%k +A ^ 


(139) 


(140) 


In the above equations, {ef} is defined in equation (37), {W} represents the corresponding 
stress resultants, Nfc is a unit global load discussed in the Section 3.3, Nya, Ny&j, and Nyajk are 
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from equation (124), and Vf n , V? t and V% are the shape functions of expanded expressions for 
node-line displacement component V n . 

Using the abbreviated notation given above, the constant coefficients of equations (134) 
can be expressed in the following form: 


p f_ 

c/= £[f {N L } T fr)dAl - J N V?\ n *dx 

p = 1 P J o n 1 

p l 

C) = V [ f ({£/&} + {Mj} T {Ki})dA-] - f N Vj'f'dx 

« J A n J 0 ‘ n. 


c iS = c> 


p ^ 

c ,y= S C f ^ N L} T G,j)dAl - J 2N y VjJ | ^dx 

P = 1 A P U 

P ^ 

c{ = V C [ 2({N*} T fty} + {/W,} r {K, y })d4] - f 2N V* | " dx 

p = 1 J * p J o * ", 

P l 

cj k - Y [ [ + {/VV} T {k,»cM] - f | "dx 


“i ^ 

c ijk = c k + cj k 


p 

<4= y [f 3{W L } r (T^}d4] - f 3W ^T' dx 

• _ J a Jn 

p= 1 * p u 1 

p ^ 

4= V [[ 3{W,}%}d/\] - f fdx 

p = 1 Ja p J ° ' "« 

p L 

Cf = V [ f 2({W, lf } r ^} + {/W^} r {K f; })dA] - [ 2 W V,? | V 

p= 1 p J ° "i 


c '"= E C J (v) r ^}d4] - 

P = 1 A P 0 i 

^ijk€ - c ijk "I" 

P L 

c} k ,= Ztf - J 

* J a Jn "i 

p= i * p u 1 


(141) 
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In the above equations, m and n 2 are the node-line numbers for the two boundary node-lines 
having edge-normal vectors in the negative and positive global y-directions, respectively. To 
perform the integration indicated in the above equations, each integrand is first expressed in 
terms of separated variables. The x-dependent terms in the integrands involve only known 
trigonometric functions, and thus are integrated analytically. The y-dependent terms in the 
integrands are evaluated at the discrete points in the y-domain of each plate strip where 
finite-difference solutions were obtained for the second-order displacement fields (see Ap- 
pendix C). One-dimensional numerical integration is then performed using Simpson's rule 
[63]. 

In the nonlinear algebraic equations (134), terms were retained through 0(q, 3 ), with terms 
of equal order in the modal imperfection amplitudes qf. This approach is helpful in determin- 
ing the form of the coefficients presented above, but it is expected that in the regime of non- 
linear response the modal amplitudes will exceed the modal imperfection amplitudes by 
perhaps one or two orders of magnitude. Assume now that 

0«,'))-0((,,) ! ) (142) 

With this assumption, equation (134) can be simplified: 

( M \ U / M \ 

xc > ~ Etf c y + + xc i ~ 

MM 1 M M M 

+ /Lj + icjjk) + X y, QflkQA c ijkf+ * c fjkA + Ofay 4 ) ( 143 ) 

;=ix = i j = ^k = ^{ = ^ 

= 0 , 1 = 1,2, ... M 

In applying equations (143), the 0(qf) contributions are neglected. 

Load-dependence of the algebraic equations: It was noted Section 3.5 that in determining 
the second-order displacement fields, a reference value must be specified for the load pa- 
rameter; the reference value is labeled Aj,. If At, is taken to be A, which is implied by some 
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researchers to be the correct approach [38], then the second-order fields and the coefficients 
of equations (141) must be computed at every load value investigated. This is undesirable 
from the standpoint of economy, and also introduces the complication that the coefficient 
equations used to determine the second-order displacement fields become singular at certain 
load values, an occurrence which is apparently a mathematical artifact not necessarily cor- 
relating to any behavior of the structure [38], Nonetheless, for the local/Euler mode- 
interaction problems investigated in the following chapter, two different methods are used; in 
the first method, X b is set to zero, whereas in the second method, A* is set to A. Results are 
compared for the two approaches. 

Solution of the nonlinear algebraic equations: As discussed in the review of literature, 
there are methods developed for following a nonlinear load-response path through limit points 
[30] and bifurcation points [31,32], and these methods provide an important capability for a 
general method of nonlinear analysis. For the test cases discussed in Chapter 4, simple 
Newton-Raphson iteration, when used with end-shortening control, has proven to be adequate 
as a solution technique. While cases were investigated which show limit-load behavior, it was 
found that limit loads are traversed in a regime in which the iterative procedure still performs 
satisfactorily using controlled end-shortening. Thus, the results presented herein were gen- 
erated using the Newton-Raphson iterative procedure. The implementation of this procedure 
is now outlined. 

For a specified set of modal imperfection amplitudes { q ° } and specified load level A, 
equation (143) can be expressed in terms of a new set of constant coefficients: 

M MM M M M 

D, + £ qPu + X X WV + X X X WtfiPw = 0 • i = 1 • 2, - M (144) 

j=1 j=V<= 1 )=u=v = 1 


where 
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( 145 ) 


M 

D,= *cf- Xcfcj 

J=J{ M 

®ij = Cjj 4 - XC ij — ^ Qk^ij 


k = 1 


D ijk = C ijk + 


XC 


ijk 


Dijkf — ^ijM 


+ XC 


A 

\}kt 


Given an approximate known solution, {q r }, representing the r* solution obtained in an itera- 
tive procedure (or the starting solution, if r = 0), define the associated residual vector, {R r }> to 
be 


M MM M M M 

R- = D, + ^ qjD (J + Qj QfPij k + y y y tf^Wfiiikt ... M (146) 


; = 1 /= 1*=1 




The residual vector of the exact solution, {0}, is expressed as a Taylor series expansion 
about the approximate, known solution: 


M 

„,i - 0 -b;+£ 


aft, 

dqj 


(qj-qf) 


{</) 


+ 0((q,-g/) 2 ) i — 1,2 M (147) 


The difference (g, — g, r ) is assumed to be small so that the square of the difference can be ig- 
nored in establishing an iterative procedure. Equation (147) can then be rearranged to get 


{g} - {q r }~ - [K r ] W 


(148) 


where the tangent stiffness matrix, [K f ], is defined by 


r dR l 
K ( = ~~z ~ 

1 dqj 


(149) 


(Q f ) 
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The Newton-Raphson iterative procedure is established using equation (148), by calling ( q } the 
(r+ 1) st solution vector: 

{q r+ ’} = { q r } - [kTW (150) 

Equation (149), when evaluated using the expression in equation (146), provides the following 
expression for the tangent stiffness matrix: 


K[j = Dlj 


M 

* = 1 


M 

fjk + Djkj) + ^ ^A^ijk + Djkjf + Dikfp 


= 1 


/,i=1,2 M 


(151) 


Convergence of the iteration procedure is assumed to occur when the following criterion is 
met: 



/ = 1 


(152) 


where e is a small positive number much less than one. 
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4.0 Results and Discussion 


In this chapter, four sets of results are presented and discussed. In the first set, the 
postbuckling response of several square unstiffened panels is Investigated. NLPAN pred- 
ictions are compared with theoretical results from the literature. In the second set, the pred- 
ictions of NLPAN are compared with test data for the postbuckling behavior of an unstiffened 
rectangular quasi-isotropic composite panel. Next, NLPAN is used to predict the local 
postbuckling response of a blade-stiffened panel of isotropic material, having three bays 
across the width. The results are compared with finite element predictions. Finally, NLPAN is 
used to predict modal interaction and imperfection sensitivity in several wide, blade stiffened 
panels, for which experimental results are available in the literature. A few comments are 
made concerning the computational cost of running NLPAN. 


4.1 Postbuckling of Square Unstiffened Panels 


4.1.1 Isotropic Plates with Simple Edge Support 


As discussed in Section 2.1.2, exact series solution procedures have been developed for 
simply supported rectangular isotropic plates subjected to uniaxial loading in the form of im- 
posed end shortening. Levy [13] obtained results for perfect square plates where the unloaded 
edges are constrained to remain straight while carrying no net edge-normal load. Coan [14] 
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obtained results for perfect and imperfect square plates for which the unloaded edges are 
uniformly free of edge-normal force resultants. In both of the mentioned works, the value used 
for Poisson's ratio was v = 0.316 . 

Predictions for end load versus center deflection computed using NLPAN are compared 
with the results of Levy and Coan in Figure 5. The end load, A/ x , is the mean value of N x over 
the width of the panel, and the center deflection is denoted w c . The eigenfunctions used in 
applying NLPAN were those deemed important in [13] and [14] (four modes for Levy's results 
and three modes for Coan's results). Thirty finite-difference intervals were used in NLPAN for 
computing the second-order displacement fields. For simple plates such as these, the method 
of NLPAN degenerates to the solution approach used by Levy and Coan. The buckling 
eigenfunctions involve only w,(x, y), and the second-order displacement fields involve only 
u, v (x,y) and v./x, y), as can be seen in equation (2). Consequently, the results of NLPAN agree 
nearly exactly with the results of both Levy and Coan, for perfect plates and for imperfect 
plates. The slight differences in results at higher load levels are most likely due to numerical 
errors associated with the finite-difference procedure used in NLPAN for determining the 
second-order displacement fields (Levy and Coan both used analytical solutions for the 
second-order displacement fields). It is interesting to note that the differences between the 
two sets of boundary conditions modelled do not affect the initial buckling problem, but do 
significantly affect the postbuckling response. 

Convergence with respect to finite-difference discretization: A convergence study was 
performed to investigate the effect of the number of finite-difference intervals, /, used In com- 
puting the second-order displacement fields. The plate reported in [13] and the imperfect plate 
reported in [14] (see Figure 5) were both modelled. The procedure used was to specify an 
end-shortening value, A u, compute the end load N x using NLPAN for a range of values of /, 
and compare N x with the values provided in [13] or [14], as appropriate. The results are pre- 
sented in Figure 6. Two postbuckling load ratios were investigated: W x /A/ xc ^2, indicated by 
open symbols, and N X IN XC ,~ 4, indicated by closed symbols, where N xc , is the critical buckling 
load. Results for Levy's boundary conditions (straight unloaded edges) are indicated with 
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Figure 5. End load varsus canter deflection for simply supported square plates: Isotropic ma- 

terial, v - 0.316 . 
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circles, and results for Coan's boundary conditions (unloaded edges free of edge-normal 
loads) are indicated with squares. 

Two observations can readily be made. First, the numerical error increases as the 
postbuckling load ratio increases. This is reasonable, since the second-order fields which 
contain the numerical error become increasingly active as the postbuckling load increases. 
The second observation is that the NLPAN results for the plate with straight unloaded edges 
exhibit much less numerical error than the NLPAN results for the plate with unloaded edges 
which are uniformly free of edge-normal force resultants. At a load of four times the critical 
value, the results for the straight-edge plate have an error of about 1.5% when ten finite- 
difference intervals are used, whereas to get the same accuracy for the other plate requires 
the use of about thirty intervals. A qualitative explanation is offered for the great difference 
in the error for the two plates. Gradients of the in-plane displacements are used in the anal- 
ysis to express the force boundary conditions at the side edges. The displacement gradients 
are expressed in terms of approximate finite-difference relationships, making them subject to 
numerical error. The plate with side edges free of load N y experiences larger in-plane dis- 
placement gradients than the plate with straight edges, and thus suffers more from numerical 
error. 


4.1.2 Composite Panel with Simple Edge Support 


Shin [64] analyzed the postbuckling response of several rectangular unstiffened com- 
posite panels with simply supported edges. The approach was essentially the same as that 
used by Levy [13] (see Section 2.1.2). For simple plate problems, the method of NLPAN de- 
generates to the solution method of [13], and thus, NLPAN should provide results in agree- 
ment with those given in [64], A single configuration in considered here. The panel is square 
with dimensions L = b = 20 inches, and h = 012 inches. The laminate configuration is 
[ ± 30l nf , where n is sufficiently large so that the laminate elastic constants Du and Du can 
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Figure 6. Convergence of results with respect to the finite-difference interval width: End 

shortening is specified, end load is computed, literature results from [13] and [14]. 
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be assumed to be zero. The lamina properties are shown in Table 2. A uniform end short- 
ening, A u, is imposed, and the y-normal edges are restrained against edge-normal displace- 
ments, i.e. v = 0. Complete details of the in-piane boundary conditions are shown In 
Figure 7. 

Table 2. Graphite/Epoxy Lamina Properties [64] 


Longitudinal Young's Modulus 

19.01 x 10 6 psi 

Transverse Young's Modulus 

1.89 x 10 6 psi 

Shear Modulus 

0.93 x 10 6 psi 

Major Poisson's Ratio 

0.38 


In both the analysis used in [64] and the analysis of NLPAN, the out-of-plane displace- 
ments can be expressed as w = X£w m „ sin(m 7 rx/L) sm{nnylb). Six modes were used in [64], 

m n 

corresponding to the coefficients Wn, w-u, Wsi, W 33 , Wji, and W 15 . Five modes were incorporated 
in the NLPAN analysis, matching those used in [64] except for the exclusion of the mode cor- 
responding to coefficient w 5 i. Results presented in [13] suggest that the excluded mode is of 
negligible importance compared to the other five modes. Twenty finite-difference intervals 
across the width were used in the NLPAN analysis for computation of the second-order dis- 
placement fields. 

The results for normalized end shortening versus normalized center deflection are shown 
in Figure 7. The critical buckling point corresponds to e, = -.00725 % (N x = -103.9 lb/in.). 
As expected, the agreement between the two analyses is essentially exact. The NLPAN anal- 
ysis predicted that the mean load N y acting on the y-normal edges is compressive for end 
shortening values up to about A u = 10 Au Cf (N x ^3AN xcr ), beyond which N y acts in tension. 


4.1.3 Square Isotropic Plates with Various Boundary Conditions 


A number of the general methods of nonlinear analysis discussed in Section 2. 2. 2. 5 apply 
only to wide stiffened panels with evenly spaced stiffeners; side boundaries are not modelled. 
NLPAN is capable of modelling a variety of boundary conditions at the side (y-normal) 
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Normalized end-shortening , au/au 



Figure 7. End shortening versus center deflection for a simply supported square unstiffened 
graphite/epoxy composite panel. 
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boundaries of a structure, as discussed in Section 3. 1.2.4 . As a demonstration of this flexi- 
bility, single mode analyses were performed for square plates using six different sets of 
boundary conditions. The exact sets of boundary conditions used are documented in 
Figure 8. Briefly described, the side-edge boundary conditions for the six sets are (listed in 
order of decreasing constraint): i) clamped, ii) simply supported (Levy's boundary conditions), 
iii) simply supported (Coan's boundary conditions), iv) one edge simply supported, one edge 
free, v) infinite width (symmetry conditions), and vi) free edges. The normalized results are 
presented in in Figure 9 and Figure 10, Results for the end load versus the center deflection 
are presented in Figure 9, and results for the end load versus the end shortening are pre- 
sented in Figure 10. While the First four sets of boundary conditions produce plates with 
postbuckling strength, the last two sets result in plates which are (approximately) neutral in 
postbuckling. 


4.2 Postbuckling of a Rectangular Unstiffened Composite 

Panel 


NLFAN was used to predict the postbuckling response of a uniaxially loaded flat unstiff- 
ened composite panel, for which experimental results are presented in [9], The test specimen 
(specimen C17 of [9]) was a sixteen ply quasi-isotropic graphite/epoxy laminate with the 
stacking sequence [±45/0/90]^. The lamina properties used in the NLPAN analysis are 
given in Table 3. The mean ply thickness given in Table 3 is based on the measured thick- 
ness of the laminate rather than the nominal ply thickness reported in [9]. The ratio of spec- 
imen length to width ( LjB ) was 3.64, and the ratio of width to thickness (B/h) was 70. The 
loaded ends of the panel were loosely clamped, and knife-edge fixtures were used to provide 
simple support slightly inboard of the two unloaded side edges. Additional details of the ex- 
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Figure 8. Labeling key for the boundary conditions considered for square isotropic piates with 
uniaxial loading. 
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Normalized compressive end load, 



Normalized center deflection, w/h 


Figure 9. End load versus center deflection for square plates having a variety of boundary con- 
ditions on the unloaded edges: isotropic material, v — 0.316 , single mode solutions 

(see key in Figure 8 for boundary conditions). 
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Figure 10. End load versus end shortening for square plates having a variety of boundary con- 
ditions on the unloaded edges: Isotropic material, v — 0.316 f single mode solutions 

(see key in Figure 8 for boundary conditions). 
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perimental configuration are shown in Figure 11(a). The test specimen was observed to 
buckle with five halfwaves over the length. Despite the fact that NLPAN cannot model the 
clamped end condition used in testing, it was felt that five was a sufficient number of buckling 
halfwaves to assure that NLPAN could accurately predict the postbuckling response at points 
away from the panel ends. 


Table 3. Graphite/Epoxy Lamina Properties [9] 


Longitudinal Young's Modulus 

131.0 x 10 5 Njcm 2 

Transverse Young's Modulus 

13.0 x 10 5 N/cm 2 

Shear Modulus 

6.4 x 10 5 Njcm* 

Major Poisson's Ratio 

0.38 

Mean Ply Thickness 

0.0125 cm (0.00492/n.) 


The details of the test boundary conditions which are shown in Figure 11(a) are not all 
described in [9], Particular omissions include the grip length over which the panel ends were 
clamped and the width of the margins between the knife-edge supports and the side edges 
of the panel. Thus, in the initial NLPAN model it was assumed that the test specimen was 
simply supported around its outer edges. The boundary conditions are given in Figure 11(b), 
and the model is referred to as Model 1. NLPAN results obtained using Model 1 showed poor 
agreement with the test data. Subsequent conversations with Marshal Rouse, co-author of [9], 
revealed additional details of the test boundary conditions which are reflected in Figure 11(a). 
A second NLPAN model, referred to as Model 2 and depicted in in Figure 11(c), was used to 
more accurately represent the test boundary conditions. Model 2 incorporates three plate 
strips, two of which are narrow strips which model the unsupported panel edges outside of the 
knife-edge supports. While Model 1 simulates the full specimen length, Model 2 simulates only 
the unsupported length of the specimen. Using Model 2, NLPAN gives greatly improved 
agreement with the experimental results (as will be presented shortly), and thus, except 
where noted otherwise, NLPAN results presented here were produced using Model 2. In the 
results presented for NLPAN Model 2, a correction is applied to the imposed end shortening 
value, which only applies to the unsupported length of the panel. The end segments of the 
panel, which are loosely clamped, are assumed to shorten as a function of the axial stiffness 
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Figure 11. Unstiffened quasMsotropic composite panel with uniaxial loading: Test description 

in (9). 
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of the unbuckled panel and the mean axial force resultant A/„, as if the stress resultant in the 
clamped ends were uniform at a value N x = N„. 

Some comments are made regarding buckling loads and mode shapes. The buckling 
eigenvalues computed using PASCO are given in Table 4 for both Models 1 and 2, and for 
halfwave numbers 3, 4, and 5. For both of the models, the ascending order of eigenvalues 
corresponds to halfwave numbers 4, 3, and 5, respectively. It is puzzling that the test specimen 
considered here exhibited five buckling halfwaves, particularly considering that the clamped 
end conditions would be expected to drive the halfwave number toward a lower value than 
what would be expected for a panel with simply supported ends. (Test specimen C18 of [9] 
was nominally identical to specimen C17, but exhibited four buckling halfwaves.) One expla- 
nation offered here is that perhaps the true buckling mode had three bulges, but in 
postbuckling an additional bulge developed at each end due to the propensity of the longi- 
tudinal wavelength to decrease with increasing load in postbuckling. 


Table 4. PASCO-Computed Buckling Loads for an Unstiffened Composite Panel 



Critical Load, N xcr ( Njcm ) 


m = 3 

3 

II 

m = 5 

NLPAN Model 1 

-880.6 

-877.4 

-961.6 

NLPAN Model 2 

CD 

CD 

x— 

« 

-1129 

-1202 


In the NLPAN predictions presented here, five buckling modes were incorporated as 
shape functions. These are the first three symmetric modes having five halfwaves over the 
length, and the first two symmetric modes having fifteen halfwaves over the length, where fif- 
teen was selected because it is three times the halfwave number of the primary buckling 
mode. (The analogous modes for a square plate correspond to the coefficients vv-n, w 13 , w 15 , 
Wji, and W 33 , where w = ][V mn sin(mirx/a) sin(n 7 ry/a), where a is the length of the sides of the 
plate. These modes were selected based on the results of Levy's analysis of a postbuckled 
square plate [13].) In computing the second-order displacement fields, twenty intervals were 
used for Model 1, whereas for Model 2, sixteen intervals were used between the knife-edge 
supports and four intervals were used for each of the outboard plate strips. In the NLPAN 
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analysis geometric shape imperfections were assumed to be present in the shape of the pri- 
mary buckling mode, at an amplitude of five percent of the panel thickness. 

Figure 12 presents the end load versus end shortening for the panel, comparing NLPAN 
predictions for both Models 1 and 2 with the test data from [9J. As can be seen, Model 2 results 
fall much closer to the test data than Model 1 results. The large difference between the results 
for the two NLPAN models shows the importance of what may appear to be minor differences 
in boundary conditions. The Model 2 results are In close agreement with the test data in the 
middle of the postbuckled range considered, but fall slightly high in the early and advanced 
postbuckling ranges. In the early postbuckling range, the test data, reported in Figure 4(d) of 
[9], may be slightly in error, since it appears that the plotted data were forced to go through 
a theoretical bifurcation buckling load which this author believes was too low. In the advanced 
postbuckling range, the test specimen has contracted in the width direction to compensate for 
large out-of-plane deflections, so that the knife-edges support the panel closer to the edges 
than in the early postbuckling regime. The effective shifting of the support toward the panel 
edges Is not accounted for in the NLPAN analysis, and may be responsible for the reduced 
extensional stiffness observed in the test relative to the NLPAN predictions. A second possible 
explanation for the disagreement in the advanced postbuckling regime is that the effects of 
transverse shear deformation, not accounted for in NLPAN, may be significant. 

Comparisons of longitudinal surface strains predicted by NLPAN with those measured 
experimentally using strain gauges are presented in Figure 13. Figure 13(a) presents the 
longitudinal strains on opposing surfaces at the center of the panel, which corresponds to a 
point of maximum out-of-plane displacement. The agreement between the NLPAN predictions 
and the test measurements is generally good. There is a small over-prediction of strain levels 
on the compression side at the higher load levels. At lower load levels, it appears that the 
agreement would be even better if a smaller imperfection amplitude were used in the NLPAN 
analysis (the amplitude used was five percent of the panel thickness). 

The distribution of longitudinal surface strains across the width of the panel at the panel 
mid-length are presented in Figure 13(b). for a fixed load level. The test data were recorded 
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End load, 



Figure 12. End load versus end shortening for a unlaxially loaded quasMsotropic composite 
panel: Test data from [9J. 
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at 3.26 times the theoretical buckling load as computed in [9], However the theoretical 
buckling load used in [9] was reported with poor accuracy (±4%), so for purposes here it 
was assumed that the theoretical buckling load in question was the value computed by PASCO 
for Model 1 for five halfwaves, which is provided in Table 4 on page 105. (This value was 
within the possible range for the theoretical buckling load given in [9], whereas the corre- 
sponding value for Model 2 was not.) A factor of 3.26 times the critical load for Model 1 (for 
m = 5) translates to a factor of 2.61 times the critical load for Model 2 (for m = 5). In view of 
the uncertainty in the load level of the test data, the NLPAN predictions are in good agree- 
ment. The subtleties of the shape of the strain distribution on the compression side of the skin 
are captured quite well by NLPAN. Note that NLPAN predicts that the highest compressive 
strains occur in the portions of the panel skin falling outboard of the knife-edge supports. 


4.3 Postbuckling of a Three Bay Blade-Stiffened Panel 


Results from the finite-element analysis of the postbuckling behavior of a uniaxially 
loaded blade-stiffened panel are presented in [4], The panel, depicted in Figure 14(a), has 
four stiffeners which delineate three bays of equal dimension across the width of the panel. 
The stiffeners are nearly three times as thick as the skin, and the panel has the approximate 
elastic properties of aluminum, given in Figure 14(a). The STAGSC-1 computer code [28,29] 
was used for the finite-element analysis; as discussed in the review of literature, this code is 
considered to be a robust method for the analysis of complicated nonlinear structural re- 
sponse. In the finite-element analysis, the panel is supported only at its longitudinal ends, 
where a clamped support condition is modelled. Uniform end-shortening is imposed. The 
early postbuckling response (through about 4 1/2 times the critical buckling load) is dominated 
by a local buckling of deflection pattern with five halfwaves in the longitudinal direction (see 
Figure 14(b)). Despite the difference in boundary conditions between the finite-element anal- 
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Axial strain, e* (%) 







ysis (clamped) and an NLPAN analysis (simply supported), it was felt that five was a sufficient 
number of halfwaves so that the finite-element results would provide a good gauge of the 
performance of NLPAN in analyzing the local postbuckling response. 

Two different NLPAN representations are used to analyze the three bay panel. The first 
is a complete (symmetric) representation, whereas the second uses a representative strip of 
panel skin centered around a single stiffener. In Section 4.3.1, results from the full panel rep- 
resentation will be discussed. In Section 4.3.2, results obtained using the two representations 
will be compared. 


4.3.1 Full (Symmetric) Representation of the Panel 


An NLPAN model was developed to analyze the response of the three bay panel. The 
response was assumed to be symmetric about the panel mid-width, so that only half of the 
panel was modelled and symmetry conditions were Imposed at the centerline. For the primary 
buckling buckling load values computed by STAGSC-1 and PASCO were N, C r = — 336.0 Ib/in. 
and N XC r = —344.4 Ib/in., respectively, where the load N x is the total end load divided by the 
panel width. (N„ here replaces of Chapter 3.) All loads reported here are normalized by 
the PASCO-computed critical load. 

In computing the second order fields in NLPAN, ten finite-difference intervals were used 
across each 4 inch skin width, and four intervals were used across each 1.5 inch stiffener 
width. A convergence study showed this level of discretization to provide solutions within one 
percent of converged values for all of the measures reported, through the highest load level 
reported. Full compatibility of displacements and equilibrium are enforced at the joints be- 
tween plate strips when computing the second-order displacement fields {u,y}. As discussed 
in Section 3.5.8, this causes the second-order fields to be dependant on the reference value 
of the load parameter, X b . It was also noted in the Section 3.5.8 that by invoking the "classical 
assumptions" of local postbuckling analysis (described In Section 2. 2. 2. 3) the second-order 
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Dimensions: Inches 

Homogeneous, isotropic material, E=l0 7 #/in^ v=0.3 
Uniaxial loading, N x = (Axial Load)/(24 in.) 

STAGSC-1 Boundary Conditions: Clamped ends with imposed end shortening 
□ •[I] - Locations of reported skin deflection and surface strains 



a) Configuration details. 



b) Planform view of primary buckling mode-shape. 


Figure 14. Three bay blade-stiffened panel configuration used for the study of local-postbuckling 
behavior: STAGSC-1 results from [4]. 
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displacement fields became load-independent. Accepting that the classical assumptions pro- 
vide accurate solutions for local postbuckling, the second-order displacement fields were as- 
sumed to be load-independent, and were computed with X b set to zero. 

Selection of buckling modes: The first issue to be taken up regards the selection of 
buckling modes for inclusion in the nonlinear analysis. Figure 15 depicts the transverse pro- 
files of the first five symmetric buckling modes having five halfwaves in the longitudinal di- 
rection. The primary buckling mode is dominated by skin buckling, with only a small amount 
of stiffener participation. Note that in the primary buckling mode the displacement amplitude 
in the center bay of the panel is about forty percent greater than the displacement amplitude 
in the outer bays. The following approach was used to select the buckling modes to be used 
in the NLPAN analysis. First, a series of NLPAN analyses was performed, each of which in- 
corporated the primary buckling mode and one of the four additional modes shown in 
Figure 15. The second and fourth modes were found to interact most strongly with the pri- 
mary mode, and were thus selected for use. Three additional modes were then selected for 
use, those being the first, second, and fourth modes having fifteen longitudinal halfwaves. 
These three additional modes have the same approximate profiles as their counterparts in 
Figure 15, except that the fourth mode having fifteen halfwaves displays much less stiffener 
participation than the fourth mode of Figure 15. The selection of secondary modes with a 
halfwave number three times that of the primary mode follows the pattern suggested by re- 
sults for the analysis of simply supported square plates [13,14]. 

Second-order field displacements at the longitudinal ends of the panel: In computing the 
second-order displacement fields for this NLPAN model in the manner described in Section 
3.5, a problem arises with regard to the satisfaction of the boundary conditions at the ends of 
the panel. The boundary conditions given in equation (29) require that at the panel ends 
(x = 0, L) the displacements v and w must be zero, neglecting the contributions of the 
unbuckled solution. For the second-order fields, however, shape functions v,y and w ti have 
contributions with an x-dependence of cos(/7mx/L), where m is an integer (possibly zero). 
Thus, v ,j and do not automatically satisfy the boundary conditions at the panel ends. The 
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characteristics of selected second-order fields computed for the three bay panel are examined 
below in order to investigate this potential discrepancy. 

The profile of the primary buckling eigenfunction (m = 5) is depicted in Figure 16, along 
with the profiles of the two contributions to its associated second-order displacement field. 
The two contributions (a = 1, 2) have values for m of 10 and 0, respectively. Scale bars show 
the approximate relative amplitudes of the two contributions to the second-order displacement 
field. For ot = 1, the profile is dominated by in-plane displacements and their gradients on each 
plate strip (a fact which is not obvious from the figure) whereas for a = 2 (where the x-de- 
pendence is constant) large out-of-plane displacements are present, clearly violating the 
simple support condition for the plate strips at the ends of the panel. The profiles of the first 
two buckling eigenfunctions having m = 5 are depicted in Figure 17, along with the profiles of 
the two contributions to the associated mixed second-order field. The comments made above 
regarding the second-order displacement field of Figure 16 also apply to the mixed second- 
order field of Figure 17. The first buckling eigenfunction for each of the halfwave numbers 
m = 5 and m — 15 are depicted in Figure 18, along with the two contributions to their associ- 
ated mixed second-order field. The two contributions have values for m of 20 and 10, respec- 
tively. In contrast to the two preceding cases, both profiles are dominated by in-plane 
displacement gradients (not completely apparent from the figure) and any violation of the 
boundary conditions at the panel ends is minimal. 

The perturbation approach used in deriving the equations governing the second-order 
fields accounted for the plate Euler equations, but the boundary conditions at the longitudinal 
ends of the structure were ignored. Obviously if the panel stiffeners were supported at their 
ends, then the displacement contributions depicted for a = 2 in Figure 16 and Figure 17 would 
be greatly altered. While a more elegant solution to this discrepancy is desirable, the simple 
fix used here to produce more physically permissible second-order fields was to impose a 
two-point constraint to all fields which are constant in the x-direction (m = 0). Figure 19 de- 
picts these constraints and their effects in modifying the profile given for oc = 2 in Figure 16. In 
Figure 19, the constraints can be seen to greatly reduce the out-of-plane displacements on 
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SECOND-ORDER FIELD PROFILES 


Figure 16. Profiles of the contributions to a second-order displacement field {<//;}. 
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BUCKLING MODE PROFILES 


\7 


x-dependence: cos(10 rcx/L ) 



a=i 


x-dependence: Constant 


a =2 




SECOND-ORDER FIELD PROFILES 


Figure 17. Profiles of the contributions to e second-order displacement field For {u,} and 

{(/;} with identical halfwave numbers (m = n). 
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Figure 18. Profile* of the contributions to • second order-displacement field {(/<>}: For {u,} and 

{uj} with different halfwave numbers (m * n). 
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each plate strip, so that the profile in Figure 19(b) is now dominated by in-plane displace- 
ments and their gradients on each plate strip. It is not apparent in Figure 19, but the in-plane 
displacement gradients on each plate strip are essentially unaffected by the addition of the 
constraints to the field. 

A few words are offered here to put into perspective the suppression of out-of-plane dis- 
placements discussed above. As was discussed in Section 2.22.2, there is justification for the 
use of the so-called classical assumptions of local postbuckling analysis [41,42], When these 
assumptions are invoked, one consequence is that w,j is identically zero on each plate strip. 
The modification to the displacement field depicted in Figure 19 accomplishes essentially the 
same thing, because the residual displacements contribute negligibly to strains compared 
to Uij and v,j. 

Results for the full (symmetric) panel representation: A comparison of NLPAN results with 
STAGSC-1 results for the three-bay panel is presented in Figure 20. For all results presented, 
the normalizing load used is the PASCO-predicted buckling load. NLPAN results are presented 
for a single mode analysis in addition to the six mode analysis. The single mode analysis was 
performed because it represents the most economical possibility for postbuckling analysis 
using NLPAN. For both load versus end shortening (Figure 20(a)) and load versus center 
deflection (Figure 20(b)), the six mode NLPAN solutions agree fairly well with the STAGSC-1 
solutions, whereas the single mode solutions exhibit a significant deviation from the 
STAGSC-1 solutions. 

The axial surface strains (e„) at the center of the panel (location 1 in Figure 14) are plotted 
as a function of the end load in Figure 20(c). The results are for a bulge at the center directed 
away from the stiffener side of the panel. In the postbuckling regime the single mode NLPAN 
solutions are skewed relative to the results of STAGSC-1, reflecting an excessive prediction 
of extensional membrane strain. The following explanation is offered. The primary buckling 
mode exhibits deflections in the center bay which are 40% greater than the deflections in the 
outer bays, whereas in the Finite element analysis the deflection amplitudes in the three bays 
were approximately equal in the early postbuckling regime [4]. The second buckling mode 
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x-dependence: Constant 



x-dependence: Constant 



Figure 19. Example of the modification used for x*independent contributions to the second-order 
displacement fields. 
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a) End toad versus end shortening. b) End load versus center deflection of skin. 



Axial strain e x (%) 

c) Axial surface strains at center of panel skin 
(location 1). 


Figure 20. 



Transverse strain e y (%) 

d) Transverse surface strains of the skin 
adjoining a stiffener (location 2). 


Analytical results for the local-postbuckling response of a three-bay blade-stiffened 
panel: STAGSC-1 results from [4]. 
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shown in Figure 15 is predicted by NLPAN to be relatively active in the postbuckling regime, 
and serves to equalize the deflections in the center bay. The excessive deflections in the 
center bay predicted using a single mode analysis induce the excessive extensional strains 
detected at the center of the panel. The six mode NLPAN solutions do not exhibit the skewing 
present in the single mode solutions, and thus give better agreement with the STAGSC-1 re- 
sults. However, the predicted strain levels are somewhat low for both the compression and 
tension sides of the skin, apparently reflecting a slightly low predicted magnitude for the 
bending curvature k*. 

Results for the transverse strains (e>.) in the skin next to a stiffener at the mid-length of the 
panel (location 2 in Figure 14) are presented in Figure 20(d). Both the single mode solutions 
and the six mode solutions show good agreement with the STAGSC-1 results up to a load of 
about twice the critical load, above which both the compressive and tensile strains are 
under-predicted. The disparity between the various predictions for transverse strains at higher 
load levels is discussed more in a following section. 


4.3.2 Panel Representation Using a Stiffener-Unit 


The three bay panel under study here was modelled using the PANDA2 computer code 
in [4]. However In the PANDA2 analysis, the panel was modeled using a repeating unit of skin 
with an attached stiffener, referred to here as a "stiffener-unit/ Since the side boundaries of 
the panel are not modelled in this representation, then an exact analysis of the finite-width 
panel is precluded, but the approach has the benefit of minimizing the complexity of the 
structural model, thus minimizing the computational expense of performing the analysis. An 
NLPAN analysis was performed using a stiffener-unit representation of the panel. This addi- 
tional analysis was performed for two reasons. The first reason was to investigate the per- 
formance of the simplified, more economical mode! in predicting the response of the three bay 
panel. The second reason was to help assess the NLPAN results obtained using the complete 
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(symmetric) panel representation. It is a straightforward matter to select the appropriate 
modes for use with the stiffener-unit model on an intuitive basis, and thus it was felt that the 
NLPAN results computed using the stiffener-unit model should serve as a check of the set of 
modes used with the complete (symmetric) representation of the panel. 

The stiffener-unit used with NLPAN is depicted in Figure 21(a). Symmetry conditions are 
imposed at each side boundary of the stiffener-unit. The load A/y is defined as the total end 
load carried by the stiffener-unit divided by the 8 inch width. While a series of three stiffener- 
units has the same width as the three bay panel, the three-unit series has less cross-sectional 
area than the panel, having one less stiffener. An area-compensated force resultant, N‘ is 
defined: 

* ZAunit (153) 

= 1.116 Wj' 

where Atotal is the cross-sectional area of the complete three bay panel, and Aumr is the 
cross-sectional area of a stiffener-unit. Therefore, in prebuckling response, equal loads N„ and 
N„* produce equal end shortening values in the three bay panel and the stiffener-unit, re- 
spectively. In presenting NLPAN results from this stiffener-unit representation, the load N x * is 
used. 

Selection of buckling modes: Five buckling modes were incorporated in the NLPAN 
analysis using the stiffener-unit representation. These were the first three unsymmetric modes 
having five halfwaves in the x-direction, the profiles of which are shown in Figure 21(b), and 
the first two unsymmetric modes having fifteen halfwaves in the x-direction. This set of modes 
is analogous to the five most important modes for an isotropic square plate, as described in 
[13]. A single mode analysis was also performed with NLPAN, because this represents the 
most economical possibility for a local postbuckling analysis. 

Results for stiffener-unlt representation: Results from one and five mode analysis for the 
stiffener-unit model, compared with the STAGSC-1 results from [4] Figure 22. Figure 22(a) 
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DEFINE: N x U =(Axial Load on Stiffener Unit)/(8 in.) 

EMU * Locations of reported skin deflection and surface strains 



i i 


STIFFENER 1 

UNIT T z 

PANEL CROSS-SECTION 

a) Stiffener unit definition. 



b) Transverse profiles of the first three unsym metric buckling modes for which m=5. 
Figure 21* Stiffener-unit representation of a three bay blade-stiffened panel. 


Results and Discussion 


123 




presents the end load versus end shortening, and the five mode NLPAN solutions can be seen 
to slightly under-predict the axial stiffness in postbuckling. Good agreement is obtained for 
the end load versus skin center deflection, shown in Figure 22(b), though this may be some- 
what of a coincidence, since there is no obvious reason why the end-load/skin-deflection re- 
sults should be better than the end-load/end-shortening results. Similarly, the fact that the 
single mode analysis seems to provide better results than the five mode analysis for end load 
versus end shortening is likely only a coincidence of cancelling errors. 

The axial surface strains in the skin at the center of the panel (location 1 of Figure 14 and 
Figure 22) are presented in Figure 23. NLPAN results obtained using the stiffener-unit rep- 
resentation are compared with those for the full-panel representation. In the single mode 
solutions, presented in Figure 23(a), the stiffener-unit model gives better results than the 
full-panel model, because the results for the former model do not exhibit the previously dis- 
cussed skewing exhibited in the results for the latter model. For the multiple mode solutions, 
presented in Figure 23(b), the two models give nearly identical results, though both results 
sets exhibit somewhat lower strain amplitudes than those predicted by STAGSC-1. 

The transverse strains in the skin surfaces at the mid-length of the panel next to a 
stiffener (location 2 in Figure 14 and Figure 21) are presented in Figure 24. NLPAN results 
obtained using the stiffener-unit representation are compared with those for the full-panel 
representation. For single mode analyses (Figure 24(a)), the full panel model provides slightly 
better agreement with the STAGSC-1 results than does the stiffener-unit model, though the 
results from both NLPAN models diverge from the STAGSC-1 results at the higher load levels. 
For the for multiple mode analyses, summarized in Figure 24(b), the two NLPAN models pro- 
vide results which agree well with each other, and which agree with the STAGSC-1 results up 
to a load level of about twice the critical value, beyond which the strain levels predicted by 
NLPAN fall below the levels predicted by STAGSC-1. 

Judging as a whole the results presented in Figure 22 through Figure 22, the stiffener- 
unit representation of the three-bay panel is fairly accurate in predicting the characterizing 
features of the postbuckling response. The good agreement between the multiple-mode strain 
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b) End toad versus cantor deflection of skin. 


Figure 22. Analytical results from NLPAN using a stiffener-unit representation of a three bay 
panel: STAGSC-1 results from [4], 
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Axial strain e* (%) 



■) Single mode solution*. 



Axial strain e* (%) 

b) Multiple mode solutions. 


Figure 23. Longitudinal surface strains in the skin at the center of a three bay panel (location 
1): Comparison of two NLPAN representations; STAGSC-1 results from [4J. 
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STAGSC-1 
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NLPAN 


Full Model - 


Stiffener Unit - 




STAGSC-1 O 

NLPAN 

Full Model 

6 Modes 
Stiffener Unit 
5 Modes 


b) Multiple mod* solution*. 

Figur* 24. Transverse aurfaca strains In the skin adjoining a stiffener for a three bay panel (lo- 
cation 2): Comparison of two NLPAN representations; STAGSC-1 results from [4], 
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predictions for the two NLPAN models (Figure 23(b) and Figure 24(b)) suggests that no im- 
portant buckling modes were omitted in applying either of the two models. It is unclear why 
the axial surface strains (Figure 23(b)) are under-predicted using both NLPAN models, even 
in the early postbuckling regime. The possibility of a reporting error in [4] can not be dis- 
counted. A question remains as to why the NLPAN predictions for the transverse surface 
strains (Figure 24) diverge from the STAGSC-1 results in the upper range of the loads inves- 
tigated. One possible explanation is offered here. The STAGSC-1 analysis predicted the onset 
of a secondary instability at a load approaching five times the critical load [4], Above this load, 
there are deflections which resemble a global postbuckling response (panel cross-sections 
translates toward the stiffener side of the panel). It is possible that even before the secondary 
instability occurs, there is load shifting occurring in the panel which is not predicted by a 
local-postbuckling analysis, and which thus modifies the deformation state compared to the 
pure local-postbuckling analysis. No global buckling modes were included in the NLPAN 
analysis because of the inability of NLPAN to model clamped end conditions. 

For the single mode NLPAN analyses performed, the stiffener-unit representation gives 
better solutions overall than the (symmetric) full panel representation, especially considering 
the strain predictions shown in Figure 23(a). As discussed before, the reason the full-panel 
representation performed poorly when a single mode was used is because the primary 
buckling eigenfunction for the complete panel is a poor representation of the early 
postbuckling displacements. A clue to this fact is that the second buckling eigenvalue is close 
to the primary eigenvalue (6% difference), indicating that in this case the second 
eigenfunction of Figure 21 plays an important role even in the early postbuckling response. 
This type of modal interaction thus may be anticipated when there is a secondary local 
buckling mode which has both the same halfwave number as the primary local buckling mode, 
and has an eigenvalue only slightly above the primary value. 
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4.4 Modal Interaction in Wide, Blade Stiffened Panels 


In this section, NLPAN is applied to wide, blade stiffened panel configurations which were 
tested experimentally in a study [25] designed to explore the interaction between the local and 
Euler buckling modes, and the resulting imperfection sensitivity. The problem of analyzing 
these phenomena in thin-walled columns and stiffened panels has received considerable at- 
tention in the literature [38,41 ,43-56,64], Several different methods have been used to obtain 
analytical predictions which are reasonably accurate when compared with test results of the 
type to be considered here from [25]. Such results are reported in the series of papers [24-27]; 
the PANDA2 panel design code was used to obtain approximate results in [65]; a method pi- 
oneered by Koiter and associates [47,52,53,55], and subsequently refined by Sridharan and 
associates [38,48,49] appears to give results in excellent agreement with test data [38]. 
Problems with perturbation approaches have been documented and discussed in the litera- 
ture [1,38,51,52], and successful approaches noted above have all abandoned the pure per- 
turbation approach. It is the intent in this section to investigate the performance of NLPAN in 
its current form (as documented in Chapter 3) which uses a perturbation approach; this hints 
that problems with accuracy may be encountered. Nonetheless, some information from the 
named references is used in this section to help guide the NLPAN analyses where flexibility 
in the method exists. In order to help illuminate some of the intuitive concepts involved in the 
analysis, the displacement shape functions to be used will be examined in detail in the fol- 
lowing section. 


4.4.1 Test Configurations 


In the experimental study of [25], several panels, each having the general configuration 
shown in Figure 25(a), were tested in uniaxial compression. Panels of two different types were 
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tested, one type having stocky stiffeners and eight bays across the width, the second type 
having thin stiffeners and nine bays across the width. The cross-sectional dimensions of the 
two panel types are shown in Figure 25(b). The panels were fabricated from epoxy resin. The 
ends of each panel were mounted in protective shoes, and the load end of the shoes featured 
a series of parallel V-grooves which spanned the width of the panel. A knife-edge supported 
each end of a panel in one of the V-grooves. The eccentricity of the load relative to the neutral 
bending axis was varied by selecting from among the parallel V-grooves in each load shoe. 
Imperfections in the shape of the local buckling modes were introduced by taking advantage 
of the temperature-dependant viscoelastic properties of the epoxy resin, A local postbuckling 
deformation shape was maintained by an applied end load while the panel was heated, then 
cooled, causing a deformation shape to set and thereby take on the role of an imperfection in 
the nominal shape of the local buckling mode. 

The load eccentricity value, e, is reported in [25] in terms of the effective Euler-mode 
imperfection amplitude, w° (the transverse deflection at the mid-length of an unloaded panel). 
In the absence of any further explanation in [25], the assumption has apparently been made 
that the two imperfection measures convert directly. In other words, it is apparently assumed 
in [25] that 

w° = e (154) 

However, the two different types of imperfections have different effects on structural response 
because an eccentric end load introduces a bending moment at the end of the panel which Is 
not present in a simply supported panel with a bowing eccentricity. An investigation of equiv- 
alency relationships for the two types of imperfections is discussed in [39], in which the fol- 
lowing relationship is recommended to relate the two measures: 

w° = 1.25 e (155) 

This issue introduces some uncertainty regarding an effective Euler-mode imperfection am- 
plitudes for the test data of (25], The local mode imperfection amplitudes reported in [25] are 
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a) Ganaral configuration of atiffanad panala. 


Dimensions: mm 


Dimensions: mm 
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based on measured deflection profiles of the panel skin along the centerline of a panel bay 
(midway between two stiffeners) of the unloaded panel. 

To model the panels using NLPAN, a stiffener-unit representation was used, as was done 
for the three bay panel in the previous section (see Figure 21(a)). Because of the large num- 
ber of bays of each panel, no correction was made to the NLPAN results to account for the 
additional stiffener present in the test panels compared to the equivalent-width representation 
using stiffener-units. The elastic properties of the panel material were not given in [28); in the 
NLPAN analysis, it was assumed that the material was linearly elastic with material properties 
E = 34.5 x 10’ Nlcm* (50 x 10 3 Ujin. 2 ) and v = 0.35 . 

The length L for each test panel was not reported directly in [25], Instead, the ratio of the 
theoretical critical loads for local and Euler buckling, PJP E , was reported. P £ is given by the 
column formula 

P E = El(j-f (156) 


where / is the moment of inertia of the panel cross-section, and P L is based in the critical 
stress for local buckling a L , given by [25] 


a L = k \ 


En 2 h 2 


12(1 — v 2 )b 2 


(157) 


where b is the bay width, h is the skin thickness, and k = 6.96 for stocky stiffener panels, and 
k = 4.0 for thin-stiffener panels. To determine the length, L, of each stocky stiffener panel for 
use in NLPAN, the equations (156) and (157) were used along with the reported ratios PJP £ 
given in [25]. Because of differences in the formulae used and the infinite-width represen- 
tation used for the NLPAN model, the ratios PJPe computed using PASCO (containing the 
VIPASA analysis) are slightly different than the nominal values reported in [25], The results for 
the stocky-stiffener panels are summarized in Table 5, along with the longitudinal halfwave 
numbers for local buckling, m L , computed by PASCO. 
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Table 5. Values of Length Computed for Three Stocky Stiffener Panels 


Panel # 

p d p E 

L 

P d P E 

m L 


(Nomina!) 

(cm) 

(PASCO) 

(PASCO) 

1 

1.50 

36.5 

1.50 

10 

2 

1.02 

30.1 

1.04 

8 

3 

0.64 

23.8 

0.68 

6 


A different approach was used to select the length for use in representing a thin-stiffener 
panel with NLPAN. The length L was simply chosen so that the ratio of critical loads predicted 
by PASCO matched the nominal value reported in [25], For the single panel considered, the 
result was L — 35.6 cm for P<./P £ = 1.05, with m L = 7. Any discrepancies between the lengths 
used in model tests and those used with NLPAN should be of little consequence in the re- 
ported results, which are presented in a normalized fashion. 


4.4.2 Displacement Shape Functions for the Thin-Stiffener Panel 


The profiles of three buckling modes of the thin-stiffener panel are depicted in 
Figure 26(a). The first mode (/ = 1) is the Euler-buckling mode (m = 1), the second mode 
(/as 2) is the primary local buckling mode, having seven longitudinal halfwaves (m = 7), and 
the third mode (/ = 3) is a secondary local mode, specifically the second unsymmetric mode 
having seven longitudinal halfwaves (m = 7). The Euler mode can be seen to involve little de- 
formation of the cross-section, and the primary local mode can be seen to involve significant 
buckling of both the stiffener and the panel skin. The secondary local mode serves a role to 
be discussed later. 

The profiles of the contributions to the second-order displacement fields associated with 
the first two buckling modes are presented in Figure 26(b). The fields depicted in 
Figure 26(b) were computed with a reference load parameter value A* = 0 (see equations 
(77)). The first pair of profiles is associated with the Euler-buckling mode ([i j] = [1,1]). Allowing 
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a) Buckling mode*. 


(i,j)=(1,1) 0-1 

[cos(2 K x/L )] 



(i.j)=(1,1) a =2 
(Constant] 


1.018 




b) Second-order displacement fields, 

Figure 26. Profiles of the displacement fields of a thin-stiffener panel: The x-dependence is 

shown in brackets. 
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for rigid-body translation in the second (a = 2) contribution, it can be seen that the out-of-plane 
displacements of the two contributions (a = 1 and a = 2) cancel for the most part at the panel 
ends, thus closely satisfying the boundary condition w = 0 for each plate strip. At the mid- 
length of the panel, the two contributions reinforce each other to produce sagging or bulging 
of the skin between the stiffeners, depending on the direction of the Euler buckling. 

The third pair of profiles in Figure 26(b) is associated with the primary local mode 
([i,j] = [2,2]). The first contribution (a = 1) is dominated by large in-plane displacements and 
their gradients, and has comparatively small out-of-plane displacements. The second con- 
tribution (a = 2), which is constant in the x-direction, exhibits extremely large out-of-plane 
displacements which grossly violate the requirement w = 0 at the longitudinal ends of each 
plate strip. As discussed in Section 4.3.1, there are grounds for suppressing these large out- 
of-plane displacements. For the results presented in following sections, constraints were im- 
posed at the node lines when computing the second displacement contribution (a = 2), in 
order to reduce the out-of-plane displacements in the same general manner as was done for 
the three bay panel (see Figure 19). 

The second pair of displacement profiles depicted in Figure 26(b) form the mixed 
second-order displacement field ([i,j] = (1,2]) associated with the interaction of the local and 
Euler buckling mode-shapes. Considering the x-dependency associated with the two profiles, 
it can be seen that the out-of-plane displacements of the two contributions tend to cancel each 
other at the panel ends (closely satisfying the boundary condition w = 0 for each plate strip) 
and reinforce each other at the panel mid-length. The net result is a deformation pattern with 
the same halfwave number as the local buckling mode (m = 7), but with an amplitude that 
varies along the length. As an illustration of this amplitude distribution, the function w 12 along 
the tip of the stiffener is shown in Figure 27. The cross-sectional profile of the mixed 
second-order field represents the modification to the local-buckling mode profile due to the 
imposition of curvature associated with the Euler buckling mode. Because the Euler mode 
curvature goes to zero at the panel ends, the amplitude of the mixed second-order field also 
goes to zero. For a positive Euler mode deflection, compression of the skin is increased and 
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compression near the stiffener tip is decreased, so the mixed second-order field serves to 
increase the out-of-plane deflections of the skin and decrease those of the stiffeners. For a 
negative Euler mode deflection, the reverse is true. Structural instability in postbuckling is due 
to a loss of bending stiffness associated with either skin postbuckling, for a positive Euler 
mode deflection, or stiffener postbuckling, for a negative Euler mode deflection. 

The second-order displacement fields {u,i} and {u K } display a negligibly small depend- 
ence on the reference value of the load parameter, 4, for the range of interest. In contrast, the 
mixed second-order field {ui 2 }, specifically its third component, w 12 , has a strong dependence 
on 4, exhibiting singular behavior at discrete values of A*. The singular behavior is depicted 
in Figure 28, which shows the maximum amplitudes of <£„ 12 (a = 1, 2) in the cross-section, as 
a function of A b . The normalizing value 4 used in Figure 28 is the critical value of end short- 
ening for the primary buckling mode. Not only is the amplitude of w 12 dependent on A b , but the 
shape is also. Figure 29 depicts the profiles of the two contributions to w i2 both for A b = 0, and 
for 4 = 0.9621. It can be seen that the bending of the stiffener section, which is pronounced 
for 4 = 0, is still present for 4 = 0.964, but is reduced in proportion to the skin displacements. 
The in-plane displacement contributions u 12 and v 12 are small relative to w 12 , and are virtually 
unaffected by the value of 4. 

4.4.3 Displacement Shape Functions for a Stocky-Stiffener Panel 


Figure 30(a) depicts the profiles of two buckling modes for the stocky-stiffener panel with 
the critical load ratio is P<./P £ = 1.02. The first mode (/ = 1) is the Euler buckling mode 
(m = 1), and the second mode (/ = 2) is the primary local buckling mode, having eight longi- 
tudinal halfwaves ( m = 8). The Euler buckling mode involves little deformation of the cross- 
section, and the local buckling mode is limited for the most part to skin buckling with little 
involvement of the stiffener. 
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Along the Stiffener Tip: w 12 = 0.08 [cos(8 jt x/L )] + 0.07 [cos(6 n x/L )] 


Figure 27. The mixed second-order displacement field at the stiffener tip for a thin-stiffener panel 

(A* - 0). 
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Maximum amplitude of <|>a 12 (Xb = ?0 
Maximum amplitude of <|) ai2 ( x 5 = 0 ) 



Figure 28, Demonstration of the singularities In the contributions to the mixed second-order 
displacement field: Thin-stiffener panel. 
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X b =0 



X b =0 96 Xi 


Figure 29. Profiles of the mixed second-order field contributions as affected by the reference 
value of the load parameter: Thin- stiffener panel. 
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a) Buckling modem. 


(UM1.1) a =1 z i : (ijHl.l) a =2 

[cos(2rcx/L )] J : [Constant] 



(ij)-(2.2) a =1 

(i.j)*(2.2) a =2 

[cos(16 Jtx/L )] 

[Constant] 

(Negligible out-of-plane 

(Negligible out-of-plane 

displacements) 

displacements) 


b) Second-order displacement fields, 3^=0. 


Figure 30. Profiles of the displacement fields of stocky- stiffener panel # 2: The x-dependence 

is shown in brackets. 
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The profiles of the contributions to the second-order displacement fields associated with 
the two buckling modes are presented in Figure 30(b), as computed with a reference load 
parameter value 4 = 0. The first pair of profiles is associated with the Euler buckling mode 
{[i,j] = [1,1]), and are qualitatively similar to the corresponding profiles for the thin-stiffener 
panel. Allowing for rigid-body translation of the second (« = 2) contribution, the out-of-plane 
displacements of the two contributions cancel for the most part at the panel ends, thus closely 
satisfying the boundary condition w = 0 for each plate strip. At the mid-length of the panel, the 
two contributions reinforce each other to produce sagging or bulging of the skin between the 
stiffeners, depending on the direction of the Euler buckling. 

The third pair of profiles, associated with the primary local buckling mode ([i,j] = [2,2]), 
are not shown In Figure 30(b), because they involve negligible out-of-plane displacements on 
each plate strip. This is explained by an inspection of the local buckling mode (/' = 2) in 
Figure 30(a). Because of the rigidity of the stocky stiffeners, each skin bay between neigh- 
boring stiffeners behaves nearly like a simple plate which is clamped at its side edges. As 
was discussed in Section 3.5.8, the second-order displacement fields contain only in-plane 
components for the postbuckling of simple plates. 

The second pair of displacement profiles depicted in Figure 30(b) form the mixed 
second-order displacement field ([),]] s [1,2]) associated with the interaction of the local and 
Euler buckling mode shapes. As was the case with the thin-stiffener panel, the out-of-plane 
displacements of the two contributions tend to cancel each other at the panel ends, and to 
reinforce each other at the panel mid-length. The complete displacement field has the same 
halfwave number as the local buckling mode ( m = 8), but the wave amplitude varies along the 
length. Here, as with the thin-stiffener panel, the cross-sectional profile of the mixed 
second-order field represents the modification to the local buckling mode profile due to the 
imposition of curvature from the Euler buckling mode. However, in contrast to the case of the 
thin-stiffener panel, the profile of the mixed second-order field here is essentially the same 
as that of the primary local buckling mode. For a positive Euler mode deflection, compression 
of the skin is increased, so the mixed second-order field iv« serves to increase the out-of- 
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plane deflections of the skin. For a negative Euler mode deflection, the compression of the skin 
is decreased, so the field w 12 serves to decrease the out-of-plane deflections of the skin. 
Structural instability in postbuckling is due to a loss of bending stiffness associated with skin 
postbuckling and a positive Euler mode deflection. A negative Euler mode deflection does not 
induce a significant postbuckling instability. Thus, imperfection sensitivity is associated with 
positive Euler mode deflections, as can be seen in the test results of [25] which are repeated 
here in a following section. 

As with the thin-stiffener panel, the second-order displacement fields {un} and {u 22 } of the 
stocky-stiffener panel display a negligibly small dependence on the reference value of the load 
parameter, X b . Component w, 2 of the mixed second-order field does exhibit a dependence on 
X b , and exhibits singular behavior at discrete values of X b , but unlike the situation with the 
thin-stiffener panel, here the cross-sectional shape of w 12 keeps its same basic shape. The 
in-plane displacement contributions u 12 and v 12 are small relative to w, 2 , and are virtually un- 
affected by the selection of X b . 


4.4.4 Analysis with NLPAN 


Sridharan and Peng [38] report analytical predictions which are in close agreement with 
the test results for the panels considered here, using a method which is highly tailored to the 
analysis of local/Euler mode interaction in isotropic columns and wide panels. Thus, while the 
method of analysis used in [38] differs In many ways from the approach of NLPAN, several 
points from the discussions in [38] are considered here for guidance in designing the NLPAN 
analyses for these panels. 

First, Sridharan and Peng state that the singularities in the mixed local/Euler second- 
order field are simply artifacts of the perturbation approach which uses modal amplitudes as 
generalized coordinates. The approach in [38] used to overcome the singularity problem is 
to include as a participating mode in the analysis a secondary local buckling mode which has 
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the same number of halfwaves as the primary local buckling mode, and which has a trans- 
verse profile similar to that of the mixed second-order displacement field. The profile of the 
mixed second-order field is then computed while constrained to be orthogonal to the second- 
ary local mode; the singularity problem is eliminated. An example of an appropriate second- 
ary local mode, consider the thin-stiffener panel, for which the profile of the secondary local 
buckling mode (/ = 3) in Figure 26(a) is similar to the profile of the mixed second-order field 
([/j] = [1,2]) shown in Figure 26(b). In the approach of [38], the thin-stiffener panel is ana- 
lyzed using all three buckling modes in Figure 26(a), which together can simulate three fun- 
damental types of deformation: i) Euler (overall) buckling deformation, ii) skin postbuckling, iii) 
stiffener postbuckling. The stocky-stiffener panel, on the other hand, would require only two 
modes: one producing Euler mode displacements, and one producing local postbuckling of 
the skin. 

Second, Sridharan and Peng use "amplitude modulation functions' (see Section 2.2 2.4) 
to allow the amplitudes of each of the local buckling modes to vary along the longitudinal co- 
ordinate of the panel. Again referring to the thin-stiffener panel, the feature of amplitude 
modulation applied to the secondary local buckling mode (/ = 3) of Figure 26(a) would create 
a displacement field which is similar in both profile and x-dependence (see Figure 27) to the 
mixed second-order displacement field. It is noted here that the phenomenon of amplitude 
modulation of a local buckling mode-shape can be simulated using multiple buckling modes 
having the same approximate transverse profile, but which differ In longitudinal halfwave 
number m. 

NLPAN analyses were performed using two different approaches. In the first approach 
(Approach 1), the second-order displacement fields were computed once with the reference 
load parameter value JLt = 0. In the second approach (Approach 2), the second-order dis- 
placement fields were computed at each load level where a solution was sought, or in other 
words, by setting = X in equation (77).. Despite the opinion expressed in [38] that the phe- 
nomenon of singularities in the computation of the second-order fields has no physical basis, 
this second approach was investigated to see if the limit point values of the panels are cor- 
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rectly predicted as the singularity is approached. All four of the panels considered were 
modelled with both approaches. The thin-stiffener panel was modelled using Approach 1, first 
with two modes, then with three, to study the importance of including the third mode. The 
thin-stiffener panel was modelled with Approach 2 using two modes. The three stocky- 
stiffener panels were modelled using two modes in both of the approaches. While it was 
mentioned earlier that the use of additional modes in the NLPAN analyses would enable the 
simulation of the amplitude-modulation phenomenon, it was felt that the important features 
of the modal interaction problem would still be captured using simple mode shapes. 


4.4.5 NLPAN Results for the Thin-Stlffener Panel 


Results for the thin-stiffener panel showing the limit load (as a fraction of the Euler 
buckling load) versus the Euler mode imperfection amplitude are presented in Figure 31 for 
two different amplitudes of imperfections in the shape of the local buckling mode. Results 
obtained from NLPAN Approach 1, both for two mode and three mode analyses, are compared 
with test results from [25], For positive Euler mode imperfections, both the two and three mode 
analysis give similar results, whereas for negative Euler mode imperfections, the two mode 
analysis does not predict any limit load below the Euler buckling load. Apparently, despite the 
presence of the stiffener buckling type of displacements in the mixed second-order field for the 
two mode analysis, the three types of displacement shapes are not sufficiently independent 
to permit the prediction of collapse induced by stiffener-dominated local buckling. Both the 
two mode and the three mode analyses under-predict the imperfection sensitivity for both 
positive and negative values of Euler mode imperfection, although the three mode analysis 
succeeds in predicting both the skin buckling and stiffener buckling modes of global collapse. 

Results for the thin-stiffener panel featuring the same test data as Figure 31 are pre- 
sented in Figure 32, but the NLPAN results presented are for the three-mode Approach 1 
analysis and the two-mode Approach 2 analysis. Both NLPAN approaches predict the two 
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Figure 31. Results for the thin-stiffener panel, NLPAN Approach 1. 
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modes of collapse (stiffener buckling and skin buckling), but both approaches generally over- 
predict the limit loads, except for large positive Euler mode imperfections, for which the test 
results show an asymptotic value of the limit load for an increasing Euler mode imperfection. 
The asymptotic trend is not exhibited by the NLPAN results. In general, the analytical results 
from Approach 2 fall closer to the test data than those from Approach 1, but considering that 
the Approach 2 analysis is at least an order of magnitude more expensive than the Approach 
1 analysis, it can hardly be concluded that Approach 2 is superior. There is irregular behavior 
of the Approach 2 results for the greater of the two local mode imperfection values, at negative 
Euler mode imperfection values (see Figure 32). This reflects numerical difficulty in computing 
the limit load when in the proximity of the singularities of the mixed second-order displace- 
ment fields. Nonetheless, it appears that both the three-mode Approach 1 analysis and the 
two-mode Approach 2 analysis are successful in representing the primary physical phenom- 
ena which produce imperfection sensitivity. 

As mentioned earlier, the local mode imperfection amplitudes reported in [25] were 
based on measurements of the panel skin deformation. However, there is no guarantee that 
the local deformations had the same transverse profile as the primary local buckling mode 
of Figure 26(a). Indeed, analytical results obtained [38] which show excellent agreement with 
the test results of [25] were obtained by taking the local mode imperfections to be in the shape 
of the secondary local buckling mode (/ = 3) shown in Figure 26(a), rather than in the shape 
of the primary local buckling mode. That choice significantly amplifies the assumed stiffener 
deformations. Making the same choice for the NLPAN analyses would result in a significant 
lowering of the computed limit loads for negative Euler mode imperfections, thus improving 
the agreement with the test data; there is, however, no information on which to base the se- 
lection of an amplitude for an imperfection in the shape of the secondary local buckling mode. 
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Figure 32. Resuits for the thin-stiffener panel comparing NLPAN Approachs 1 and 2. 
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4.4.6 NLPAN Results for the Stocky-Stiffener Panels 


Results for the three stocky-stiffener panels (see Table 5 on page 133) are presented in 
Figure 33 through Figure 35. Results from two mode analyses, using both Approach 1 and 
Approach 2, are compared with the test data from [25]. The panels were each tested at two 
different local-mode imperfection amplitudes (except for Panel 3 which was tested at a single 
local-mode imperfection amplitude) over a range of Euler-mode imperfection amplitudes. 
Results from the two NLPAN approaches agree fairly well with each other, except for Panel 3 
for negative values of imperfections in the Euler mode-shape. In this latter case, the irregular 
results of Approach 2 are apparently due to numerical problems associated with the proximity 
of singularities in the mixed second-order Field. For all three panels, there is fair qualitative 
agreement between NLPAN predictions and the test data for small amplitudes of the Euler 
mode imperfection, but the test data for all three panels exhibit an asymptotic value of the limit 
load for an increasing amplitude of the Euler mode imperfections, and NLPAN fails to predict 
this phenomenon. The results for the two NLPAN approaches match each other closely 
enough that the extreme economy of Approach 1 relative to Approach 2 makes the former 
approach the preferred one. 


4.4.7 Additional Comments on the NLPAN Results for Imperfection 
Sensitive Panels 


There are several factors which bring uncertainty into the attempts to correlate the ana- 
lytical predictions of NLPAN with the test data of [25]. These factors include the use of linearly 
elastic material properties to model the epoxy material, the use of bowing Imperfections to 
model load eccentricities, and the use of a single local buckling mode-shape to model the lo- 
cal shape imperfections. In a related concern, Sridharan and Peng [38] reported that in order 
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Figure 33. 


Results for stocky- stiffener panel # 1 comparing NLPAN Approachs 1 and 2. 
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Figure 34. Resuits for stocky-stiffener panel # 2 comparing NLPAN Approachs 1 and 2. 
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Figure 35. Results for stocky-stiffener pane! # 3 comparing NLPAN Approachs 1 and 2. 
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to get good agreement with the test data for the stocky-stiffener panel with the critical load 
ratio PJPe = 0.64, it was necessary to incorporate an expression for the mid-surface curvature 
of plate strips, k x , which accounts for large rotations associated with the large Euler mode 
deflections which are encountered. While better accounting for some of these factors might 
improve the agreement between the present analysis and the experiments, there seems to 
be a fundamental shortcoming in the NtPAN approach so that the asymptotic behavior of the 
limit load with respect to increasing Euler mode imperfections is not predicted. (The presence 
of a lower bound to the limit load observed in the test data for positive Euler mode deflections 
was also predicted analytically for stocky-stiffener panels by Koiter and Pignataro [52].) 

Of the two NLPAN approaches used, Approach 1 is judged to be superior for reasons of 
economy. However, as seen with the thin-stifTener panel, it is important to select the proper 
combination of local buckling modes to assure that all dominant type of deformation can vary 
independently. A theoretical drawback of Approach 2, not discussed earlier, is that the pres- 
ence of singularities in the the second-order displacement fields at certain critical values of 
X produce singular behavior in the predicted structural response as the end shortening X 
traverses those critical values. Thus, while Approach 2 produces limit load predictions com- 
parable to those determined from Approach 1, the prediction of continuous structural behavior 
as the end shortening increases beyond the load-limit point is precluded in Approach 2. 

The various methods which provide reasonably accurate results for these imperfection 
sensitive panels (described in the references listed at the beginning of Section 4.4 all make 
use of analytical models which place various restrictions on the configuration geometry. These 
restrictions prohibit the direct incorporation of any of the methods into the general method of 
NLPAN, because one of the goals here is to retain the configuration flexibility of the buckling 
analysis in V1PASA. Improved accuracy of the NLPAN method in analyzing local/Euler mode 
interaction is desirable. Should further attempts to refine the theory behind NLPAN fail to 
produce the desired improvements in accuracy, then it would be worth considering imple- 
menting concepts from some of these other methods into NLPAN for application to the limited 
classes of configurations to which they apply. 
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4.5 Computational Expense of NLPAN 


The version of NLPAN used to generate the reported results is still in a state of develop- 
ment, and thus, a number of refinements which would improve the computational efficiency 
have not yet been implemented. Nonetheless, some execution times are given here in order 
to characterize the computational expense of NLPAN. The times reported are for runs made 
on an IBM 3090-300 digital computer. 

Computational costs will be discussed for the symmetric representation of the three bay 
blade stiffened panel (see Section 4.3.1), because this model has the most complicated 
cross-section of any of the configurations investigated. The model features four plate strips 
and five node lines. A total of 42 discrete y-coordinates were used on the cross-sections of 
the various plate strips, both for obtaining the finite-difference solution to equations (96), and 
for performing numerical integration in the y-direction in evaluating the area integrals ap- 
pearing in equations (141). The central processor unit (CPU) times required by NLPAN to 
perform 1, 2, 3, 4, 5, and 6 mode analyses were 1.7, 2.8, 5.5, 12, 25, and 48 seconds, respec- 
tively, not including the run times for the PASCO (VIPASA) analyses (used to provide the 
buckling mode shapes), which were small in comparison. Note that the NLPAN run time is 
approximately doubled each time an addition mode shape is incorporated in the analysis. 

To get an indication of the improvements in computational economy which can be real- 
ized by implementing the "classical assumptions' for local postbuckling analysis (see Section 
2.2. 2.3), an additional NLPAN run was made with the following modification: the integrals ap- 
pearing in equations (141) were simplified by incorporating the assumption that 
u l = vi = w lj = 0 (/,/ = 1, 2, ... , M). (The indicated displacement contributions are zero as an 
outcome of imposing the classical assumptions.) For a 6 mode analysis, the modification de- 
creased the CPU time from 48 seconds to 17 seconds. The results of the analysis were es- 
sentially unaffected by the modification, supporting the validity of the classical assumptions. 
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5.0 Conclusions and Recommendations 


5.1 Conclusions 


A method has been presented for the geometrically nonlinear analysis of compressively 
loaded prismatic composite structures. The permitted structural loading is uniaxial com- 
pression, or for flat panels, biaxial compression. Geometric shape imperfections are permit- 
ted, and the longitudinal ends of the structure are assumed to be simply supported (transverse 
deflections are zero). Structures are modelled as assemblages of flat plate strips which are 
rigidly joined along mutual longitudinal edges. In order to permit the accurate modelling of 
various structural sections, small eccentricities can be simulated between the reference edge 
lines of adjoining plate strips. 

The method is semi-analytical in nature, and thus, is computationally economical when 
compared to finite element methods. Buckling eigenfunctions, determined using the VIPASA 
computer code (the primary analysis code within the PASCO sizing code for stiffened panels), 
are used as the primary shape functions for the displacements in the nonlinear analysis. The 
modal amplitudes associated with the buckling eigenfunctions serve as the generalized co- 
ordinates in a nonlinear analysis. Displacement contributions which are of second order in 
the modal amplitudes are also incorporated. The principle of virtual work is used to obtain 
an approximate expression for the equilibrium condition, which takes the form of a set of 
nonlinear algebraic equations involving a finite basis of modal amplitudes. 

The method has been implemented on a digital computer in a FORTRAN computer code 
designated NLPAN. Central processor unit times for the NLPAN runs were found to double 
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with each addition mode shape incorporated in an analysis. For situations where they are 
applicable, implementing the "classical assumptions' of local postbuckling analysis into 
NLPAN would provide large reductions in the computational cost of the method. 

While the method is applicable to a wide variety of structural configurations, those in- 
vestigated here were selected based on the need during the development of the method to 
start with simple configurations and systematically increase the complexity. The verification 
test cases include several stiffened and unstiffened panel configurations for which results are 
available in the literature. NLPAN was demonstrated to provide accurate results for the local 
postbuckling behavior of unstiffened and stiffened panels through moderately large 
postbuckling loads. 

For blade-stiffened panels exhibiting interaction between local and global buckling 
modes, NLPAN is successful in predicting sensitivity to geometric imperfections in the shape 
of the local and/or global buckling modes. However, both theoretical and experimental results 
in the literature demonstrate that for certain types of stiffened panels the loss of strength as- 
sociated with a global-mode imperfection has a limiting value for an increasing imperfection 
amplitude, and this phenomenon is not predicted by NLPAN. This apparently reflects an in- 
herent limitation of the perturbation approach used in the current method. 

Two different weighted-residual approaches were explored for obtaining a set of nonlin- 
ear algebraic equations governing equilibrium. However, it was discovered that there are 
unbalanced force- and moment-resultants at the edges of the plate strips and along joint lines, 
which were not accounted for in the weighted-residual approaches, thus resulting in gross 
violations of the virtual work statement (for some configurations). It was decided to employ the 
virtual work statement directly for obtaining the nonlinear algebraic equations. 

The assumed functional form for the second-order displacement fields was chosen so as 
to permit separation of the variables x and y in the differential equations governing the fields. 
This has the consequence that the boundary condition at the longitudinal ends which requires 
that the transverse deflections be zero is not automatically satisfied. In applications to the 
verification test cases, however, it was found that the second-order displacement fields do 
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closely satisfy the stated boundary condition, except for contributions which are both associ- 
ated with local postbuckling and are x-independent. The difficulty with this one special case 
could be eliminated with the adoption of the so-called "classical assumptions" for local 
postbuckling analysis (see Section 2. 2. 2. 3), the use of which appears to be well justified for the 
configurations explored. 

The equations governing the second-order displacement fields, determined using a per- 
turbation expansion of the plate equilibrium equations, possess a term containing the end- 
shortening control parameter, suggesting that the second-order displacement fields should 
be load dependent. However, that is highly undesirable from the standpoint of computational 
economy. For the analysis of local postbuckling, the use of the "classical assumptions" men- 
tioned above removes the load-dependence of the second-order fields, and thus it was as- 
sumed here that for local postbuckling analysis the second-order displacement fields could 
be considered to be load-independent. For the analysis of local/global model interaction, the 
second-order fields associated with global buckling were found here to be only slightly load- 
dependent over the load range of interest. On the other hand, the mixed local/global second- 
order fields are highly load-dependent, both in amplitude and in shape, and even become 
unbounded in amplitude for certain values of the load parameter. This singular behavior 
poses a fundamental problem with the use of load-dependent second-order displacement 
fields. As reported in the literature, this difficulty has been dealt with and overcome in ap- 
proaches which consider some special classes of configurations, but whether analogous fixes 
could be made to the general method presented here is a question which requires further 
study. The most successful approach attempted here for the analysis of local/global mode 
interaction was to compute the mixed second-order field at a zero load level, and add as a 
participating mode a second local buckling mode with the same halfwave number as the pri- 
mary local mode, but with a transverse profile similar to that of the mixed second-order dis- 
placement field. As mentioned previously, there remain accuracy problems with this 
approach. 
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Results for locally postbuckled rectangular plates suggest that the primary buckling mode 
is a good indication of the deflection shape in the early postbuckling regime. Furthermore, the 
selection of additional modes to be included in an analysis is based on the idea of refining the 
dominant mode shape while retaining its basic symmetries or anti-symmetries. However the 
selection of modes for use in the analysis of more complicated structural sections was found 
here to be less straightforward. In a local postbuckling analysis of a blade-stiffened panel with 
three bays across the width, it was found that although the primary buckling mode appeared 
to have the general shape expected for the early postbuckling displacements, a second 
buckling mode, having the same longitudinal halfwave number as the primary buckling mode, 
played a significant role in representing the displacements of early postbuckling. The sec- 
ondary buckling mode mentioned had a critical load only 6% greater than the primary 
buckling eigenvalue, suggesting that a close proximity of eigenvalues should be taken as a 
hint of this particular type of modal interaction. 


5.2 Recommendations for Future Work 


Several recommendations are made here for improvements to and additional assess- 
ment of the method of NLPAN. These are: 

1. Incorporate the so-called "classical assumptions" of local postbuckling analysis, thereby 
improving the computational economy and eliminating certain discrepancies in the satis- 
faction of boundary conditions at the longitudinal ends of the structure. 

2. Explore modifications to the theoretical approach which will improve the accuracy of the 
method in predicting the response of structures exhibiting local/global mode interaction. 
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3. Investigate the performance of NLPAN in analyzing the response of configurations which 
are subjected to biaxial loading. 

4. Investigate the performance of NLPAN in analyzing the response of shell-like structures, 
specifically structures featuring plates which join at shallow angles. 

5. To improve the computational economy, explore the practicality of replacing certain nu- 
merically derived solutions used in NLPAN with analytically-derived solutions. The par- 
ticular solutions of concern are the solutions for the transverse-dependency of the 
second-order displacement fields, and the solutions for integrals which are evaluated 
over the cross-sectional domain of the structure in determining the coefficients of the 
nonlinear algebraic equations. 

6. Implement continuation methods into the solution procedure to permit the negotiation of 
limit points and secondary bifurcation points. 
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Appendix A. Nonlinear Plate Theory 


The plate theory developed here applies independently to each plate of the linked-plate 
structure. While large displacements may occur, material strains are assumed to be small (as 
limited by material failure) and Hooke's law is assumed to govern material response. Plates 
have elastic properties consistent with the assumptions of classical laminated composite plate 
theory [61]. The Kirchhoff-Love assumptions are used for specifying the distribution of dis- 
placements through the thickness of the plate. Geometrically nonlinear effects are accounted 
for in a manner consistent with the plate theory of von Karman [66], but with certain modifi- 
cations required due to the nature of displacement fields experienced by certain component 
plates in a linked-plate structure. The principle of virtual work is used to derive the boundary 
value problem governing equilibrium of the plate. A development of the plate theory follows, 
which includes a consideration of geometric imperfections. 


A.1.1 Strain-Displacement Relations 

Let {u} = [u(x, y) v(x, y) w(x, y)] r denote the three displacement components of the plate 
mid-surface in the x-, y-, and z-directions, respectively, where the x-y plane lies at the mid- 
surface of the undeformed plate (see Figure 36. Fora Lagrangian description of deformation, 
the finite-strain expressions for the in-plane strain components of the mid-surface are given 
by [67] 
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H = U ’X + y ( U 'X + V 'X + w ’x) 

Ey = V, y + (U,y + V,y + VV,y) 

y^ y = U.y + V, X + U, X U,y + V, X V,y + W, X W,y 

In the plate theory of von Karman, only displacement gradients w„ and w, y are expected to 
achieve significantly large amplitudes, so of the nonlinear terms in equations (A1), only 
w,| , w, £ , and w,„w,j, are retained. However, in the present application, gradients of u and v, in 
addition to gradients of w, may become significantly large due to in-plane rotation of a com- 
ponent plate, as would occur in a stiffener whlcTT participates in the global buckling motion of 
a stiffened panel. Thus, displacement gradients due to in-plane rotations are explored here. 

A 

Consider a local rotation with a vector denotation of co = cofc using the right hand rule, 
where ^ }s the rotation amplitude, and k is the unit vector in the z-direction. The in-plane 
displacement gradients corresponding to this rotation are expressed in terms of Taylor series 
expansions in co: 

= v <y = - (1 - cos to) = - -y- + O(t 0 4 ) 

— u ty = v, x = sin co = co + 0(o> ) 

Displacement gradients u }X and v, y are of order co 2 , whereas u, y and v tX are of order co. It is 
assumed that co is sufficiently small to permit the neglection of terms of order (a> 3 ) in equations 
(A1). The following expressions for the mid-surface In-plane strain components are obtained: 


£ x = u 'x + \ ( v 'l + w «x) 

Ey = V,y 4 - — ( U,y + W.y) 

y°xy “ u »y + v *x + w 'X w *y 


In the above equations, the terms v,? and u} y represent contributions not found in the von 
Karman nonlinear strain expressions. 
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The Kirchhoff-Love assumptions are used to determine the distribution of the in-plane 
strain components through the thickness of the plate. The assumptions are that material lines 
perpendicular to the mid-surface of the undeformed plate remain straight, inextensional, and 
perpendicular to the mid-surface. Donnell [66] derives the following relations for the distrib- 
ution of strains through the thickness, using the Kirchhoff-Love assumptions: 

e x = U, x + y (v,l + w 2 x ) + z[ - w, xx + u, xx w, x + v, xx w, y + u, x w , xx : ] 

By = V,y + y (U 2 + W 2 ) + Z[ - W, yy + U,yyW, x + V,yyW,y + V. y W, yy ] (A 4) 

Yxy = ~ u ~x V ’X ~ u <y v <y 

+ Z[ - 2W, Xy + 2U, Xy W, X + 2V, X yW,y + 2 (U, X + V,y)W, X y + ( U,y + V ' X )(w ' XX + W, yy )] 

where terms of higher order and terms nonlinear in z have been neglected. If terms in 
equations (A4) which are independent of z (the mid-surface strains) are retained or discarded 
using the same assumptions as used in deriving equations (A3), then the mid-surface strains 
of equations (A4) reduce to equations (A3). It is further assumed that the plate thickness and 
the displacement derivatives are all sufficiently small so that only the lowest order z-depend- 
ent contributions need to be retained. The resulting z-dependence is the same as that of 
classical linear plate theory. The final expression for the distribution of in-plane strains in the 
plate is given by 

e x (x, y, z) = e x (x, y) + zk x (x, y) 

e y (x, y, z) = e y (x, y) + zx y (x, y) (A 5) 

Y X y(x, y, 2) = V C xy ( x . y) + ZK xy( X ' Y) 

where the mid-surface strains are given in equations (A3), and where the mid-surface curva- 
tures are given by 

K X = — W 'XX 

Ky= - W iyy (A6) 

K X y ~ — 2 W, X y 
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In the presence of geometric shape imperfections, it is assumed that when the imperfect 
plate is unloaded and free of stress resultants, the displacement variables describe the 
imperfection shape: {(/} = {u 0 } = [u° v° w°] r . In order to establish the "mechanical strain," 
meaning the change in strain due to the application of structural loading, express the mid- 
surface strains and curvatures, {e c } = [ej t c yy%] T and {k c } = [kS k c y Kj y ] r , respectively, as func- 
tions of the displacement components: 


{e c } = {e c ({u})} 
{k c } = {k c (w)} 


m 


where the functional relations are given in equations (A3) and (A6). The mid-surface me- 
chanical strains and curvatures, {e m } and respectively, are given by the difference be- 
tween the apparent strains and curvatures (those based on {u}) and the initial strains and 
curvatures (those based on {u 0 }): 

rlu, x + .5(v, 2 x + w*)] - [u° + ,5(v,° 2 + w? x )]-j 
(e m ) = {e c ({u})} - {e c ((u 0 })} = < [v, y + .5 (u, y + w, y )] - [v° + .5 {u?* + )] l (A8) 

Uu.y + V, X + W, X W,y] - lu°y + V? X + W° X W?y ] J 

CW, xx -w? xx -j 

{k" 5 } = (K C (W)} “ {-cV 0 )} = - S ^,yy - W?yy > (A9) 

(-2W, xy - 2W° y J 

The other three stress/strain components are considered briefly for later use. The 
Kirchhoff-Love assumptions imply zero transverse shear strain, and this condition is ex- 
pressed as 

Yxz = Yyz = ° 


Transverse normal stress is assumed to be negligible compared to in-plane stresses, so the 
assumption is made that 

<r z = 0 (AH) 
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where o z Is the normal stress perpendicular to the plate mid-surface. 


A.1.2 Equilibrium Equations in Terms of Stress Resultants: 

Consider a linearly elastic body In the absence of body forces, undergoing large dis- 
placements but small strains. The principle of virtual work (equivalent, for this conservative 
system, to the principle of minimum total potential energy) provides the following equilibrium 
condition [68], where tensor notation Is used: 


f a,jS tijdV- f t/Su/ dS-0 i, j = 1, 2, 3 (A 12) 

**V 

where summation on / and j is implied, and where S 2 is the portion of the boundary surface 
on which tractions, t, are specified, Sdj is any kinematically admissible virtual strain field 
within the domain of the body, and Su t is any kinematically admissible virtual displacement 
field for the boundary surface. On portions of the boundary surface where tractions are not 
specified, displacements are assumed to be specified. 

For application to plates, integration over the volume is replaced by successive inte- 
grations over the plate area, A (in the x-y plane), and through the plate thickness, h (in the 
z-direction). Integration over the boundary surface is replaced by successive integrations 
along the tangential edge coordinate, s, and through the plate thickness, h, where it is as- 
sumed here that all non-zero boundary tractions are applied at the plate edges. The stress 
and strain tensor components are converted to the symbols preferred here: 

[<*11 • <*22 < <*33 ■ <*23 • <*13 • <*12^ = [<*x • <*y ■ <*Z > T yz • T xz • T xy^ 

(413) 

C«ii ' e 22 > ^3 ■ ^ e 23 . 2e^3 , 2e^] = > E y ■ e z < fyz < Txz • Txy] 
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where it is recognized that the stress and strain tensors are symmetric. Accounting for the 
zero stress and strain components specified in equations (A10) and (All), the virtual work 
expression can then be written 

Jl Jl 

n 2 (o x Sc x + <7 ySty + T xy Sy xy ) dzdA- j | t i 5u i dzds = 0 i = 1, 2, 3 {A 14) 

2 2 

where S 2 is the portion of the edge coordinate along which tractions are specified on the edge 
surface. 

Now define generalized force resultants on the boundary, f x , f y , f z , and M„, which are de- 
picted in Figure 37(a). Resultants fj y , and f z are edge forces per unit tangential length di- 
rected parallel to the coordinate axes. Resultant M„ is the edge bending-moment per unit 
tangential length, and its line of action follows the plate edge during deformation, so that the 
vector direction of the moment is not constant. The edge twisting moment M , , shown in 
Figure 37(b), is not included in the above group of stress resultants, because f z is assumed to 
contain the Kirchhoff equivalent shear force due to the edge twisting moment. 

The previous definition of S 2 is refined here. Assume that prescribed boundary tractions 

* A A A 

are specified in terms of the four generalized force-resultant values f Xl f yi f Zt and M ny each of 
which acts on a specific portion of the boundary. Denote these four different portions as 
Ci,C 2 , C 3 , and C 4| where these portions may overlap or coincide. When the strain forms of 
equations (A5) are substituted into the virtual work expression (equation (A14)) and the inte- 
gral through the thickness is evaluated, the resulting expression is 


f (N x 6t C X + NySZy -h N X ySy x y + ^X^ K X'XX + MydKy t yy + M X y 6 K X y, X y) dA 
J J 


- j f x Su ds - f y Sv ds - f z Sw ds + M n Sw, n ds 
in. Jc* 


= 0 


where the stress resultants are defined by 


(A 15) 
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a) Generalized edge force-resultant# used in the nonlinear plate theory. 



Figure 37. Measures of the loads acting on a plate edge. 
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The stress resultants are depicted in Figure 36. 

Substituting the expressions for mid-surface strains and curvatures of equations (A3) and 
(A6) into equation (A15) and applying the two-dimensional form of Green's Theorem, the area 
integral of equation (A15) can be converted into the following expression: 


— f {CNx'X 4" N X y,y + ( NyU,y),y]SU + (^X^'x)«X^^ 

J A 

+ + 2M X y, X y + My.yy + (N^,* + N X yW,y), X + (NjfyW,* + Ny W , y ) , y ] 6 W} tf A 

+ I* {t n x N x + n y( N xy + W y u, x )]<3u + ln x (N X y + N x v, x ) + r>yNy]5v 

J s 

+ ln x (M X , X + M X y,y + N X W, X + N X yW,y) + + My,y + N Xy W, x + N y W, y )]<5W 

— (n x M x + riyM X y)6w, x — (n x M xy + n y M y )Sw, y } ds 

where n x and n y are the components of the unit normal vector of the undeformed plate edge, 
defined by 


A A A 

n = n x i + n y j 


( 418 ) 


It is convenient at this point to refer certain terms to a natural edge coordinate system 
(ns), where n and s are directed along unit vectors n and s, respectively. Vector n is defined 
in equation (A18) and s is defined by 

s = -n y i+n x j (419) 

These unit vectors are depicted in Figure 37(a). Derivatives in x and y can be expressed in 
terms of the natural edge coordinates: 


Appendix A. Nonlinear Plate Theory 


172 



(A 20) 


_j d_ n _5_ n __d 
dx x dn y ds 

-A-n 5 I n _ A 

dy y dn x ds 


Using these relationships, the last three terms in expression (A17) can be converted as fol- 


lows: 


{[n x (M x , x 4 - M xyiy + N x w, x 4 - N xy w,y) 4 - n y (M xy , x 4 - M yiy 4 - A/ xy w, x + N y w ty )]Sw 

- (n x /W x + n y M xy )5w, x - (n x M xy + nyM y )8w, y } ds (A21) 

= I {(Q n + A/„w, n + N s w, s )<5w - M n Sw tn - M s <5w, $ } ds 
J s 

where the newly introduced stress-resultants are given by the following expressions [68]: 

Q n = (M X , X + M xyt y)n x + (^Xy’X + My }y )n y 

N n = N x”x + N y n y 4- 2n x nyN xy 

= (N y N x )n x n y ^xy(^x ^y) (A 22) 

= M x n x 4- M/i y 4* 2n x ny/W xy 
M s = (Af y - M x )n x n y + M xy (n x - nj) 

The stress resultants defined above are depicted in Figure 37(b). The last term in equation 
(A22) is now manipulated to establish the Kirchhoff equivalent shear force contribution of M s . 
If we require that the edge-tangent vector s varies continuously with s except at coordinates 
Si, s 2l ... s N , then applying Green's Theorem it can be shown that _ 


J * f /I s ! I *2 I 5 /V 

— M $ 6w lS ds = M s , s 6w ds 4- M s 6w _ -f M s Sw I + ... 4- M $ Sw 
s V *1 *2 S N 


The terms in parentheses represent the virtual work due to concentrated forces at Si, s 2 , ... s w , 
which are usually neglected in applications of plate theory, as they will be here. 
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The equilibrium condition is now reassembled using the substitutions developed above. 
It is noted that for edge coordinates not lying on the respective boundary portions 
Ci, C 2 , C 3 , and C 4 , the corresponding virtual displacements Su, Sv, Sw, and Sw,„ , are zero. 
Thus, the virtual work statement (equation (A15)) can be written as 


— | (EMjc’X "b ^xyy "b (NyU<y),y~]$U + EWxyx "b ^yy "b (^x^’x)’x^^ 

*b x 'xx + 2Mxy*xy + Afy,yy "I - (A/^. w,^ 4- N xy W ty ) , x “b {N xy w, x H“ N y w , y } , y ~\Sw] c 1A 

+ ln x N x + n y (N xy + N y u, y ) - f x lSu ds 4- ln x (N xy + N x v, x ) + n/J y - f y ~\6v ds 

J c , J Cj 

+ | [Q n + M s , s + N n w, n + N s w, s - f z ](5w ds - [ lM n Sw, n - M n ]Sw, n ds 
Jr. J C. 


(A 24) 


= 0 


By the fundamental lemma of the calculus of variations, Su, Sv, and Sw may vary inde- 
pendently, and may vary within the domain independently of the values at the boundaries. This 
supplies the Euler equations by requiring that the bracketed expressions in the area integrals 
of equation (A24) be identically zero over the domain of the plate: 

^X’X "b A/ X y,y + (NyU, y ), y = 0 

tycyx "b ^yy "b (^x V 'x')’x = ® (A 25) 

M X’XX "b 2A ^xyxy *b ^ yyy "b (N X W, X + A/ X yW,y), x 4- (W X yW, x 4- A/yW,y),y — 0 


The boundary conditions are obtained by requiring that each boundary integral of equation 
(A24) be zero. One possible set of boundary conditions is obtained by requiring that the fol- 
lowing force or displacement boundary conditions must be satisfied at all points of the 
boundary (plate edge): 


A 

h = 

or 

c 

II 

c:> 

II 

0 s 

or 

A 

V — V 

<0< 

II 

O' 

or 

<i t 
II 
* 


M n = M„ or w,„ = w, n 


(A 26) 
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where expressions for edge force-resultants /*, % . and f z are extracted from equations (A24): 


f X — n x^X n y(^xy + NyU,y) 

fy = n x (N xy + N x v, K ) + n y N y (A27) 

f z = Q n + M s , s + N n w, n + N s w, s 

Edge force-resultants f x and f y are not expressed in terms of the stress resultants defined in 
equations (A22) because of the complicated, counter-intuitive form which results. 

In the current application, all plate strips are rectangular in shape, so that all plate edges 
can either be classified as x-normal (n = ± i) ory-normal (n = ±j). Equations (A26) and (A22) 
can be evaluated to obtain, for an x-normal edge. 


K — n x N x 

fy ~ n x(^xy ^x V ’x) 

f z = n x {M X , X + 2M X y,y + N X W, X + N X yW,y) 

M n = M x 


(A 28) 


and for a y-normal edge. 


f X — riy(N X y + NyU,y) 
fy " ' r?yA/y 

f z = n y (2M X y, X + My,y + N X yW, X + NyW.y) 
M n = My 


(A 29) 


A.1.3 Plate Constitutive Equations: 


In the conventional notation of laminated composite plate theory [61] the stress resultants 
and the mechanical strains and curvatures at the mid-surface have the following relationship: 
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( 430 ) 
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where the square matrix is the stiffness matrix of the laminate. All elements of the stiffness 
matrix are constant for a given laminate, and follow conventional definitions given in, for ex- 
ample, reference [ 61 ]. 
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Appendix B. Y-Dependence of Buckling 
Eigensolutions 


In this appendix, the functional forms for the buckling eigenfunctions are developed. The 
basis for the development presented here is the development provided in [5] pertaining to the 
VIPASA computer code. Much of the notation used here follows that of [5], though some 
changes were made to accommodate notation preferred in the current development. Aspects 
of the development in [5] which are related to vibration or to loading and plate anisotropic 
constants not being treated by the current method, are omitted here. The focus in [5] was on 
obtaining stiffness matrices for plate strips, so the functional forms for buckling eigensolutions 
were not provided. Thus, the functional forms are derived here. 


B.1.1 Complex Representation 


Variables x and y describe the mid-surface coordinates of a plate strip, where x lies in the 
interval [0, L], and y lies in the interval [0, b]. The integer m denotes the number of halfwaves 
over the length L corresponding to the buckling eigenfunction under consideration, and 
dimensionless variables X and Y are introduced, defined by the expressions 

X -HLEJL Y = -^~ (B1) 
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A buckling eigenfunction, here denoted as [u vw], is represented using complex functions of 
X and 7. The following form is assumed in VIPASA, for the special case when the plate elastic 
constants D 16 and D 26 are zero: 


u(x,yK 
v(x, y) t = Re 
w(x, y)J 


U(Y), 
V(Y) L ,x 
W(Y)J 


(62) 


where /' = J-\ , the functions U, V, and W are in general complex, and Re signifies the real 
portion of the complex expression in the brackets. For the current application, a phase shift 
Is applied to the assumed form of equation (B2), by taking as the displacements the imaginary, 
rather than real, portions of the complex functions. This is done to locate the x-domain of the 
structure in the interval [0, L]. The new expression is 


u(x,yK 
Kx, y) > = Im 
w(x, y)J 


U(Y), 
V(V) >e ,x 

W(y)J 


(63) 


The complex functions [L/ V W] are represented in terms of two sets of real functions, 
lu„ V R W„] and [U, V, W,]. as shown: 


ri'OO', rU R (Y), rU,(Y), 

^(y) U } v r (Y) U i) y,(y) l (B4) 

lw(y)J w R (y)J ^vv^voJ 


In the current application the applied shear loading and the elastic constants Die and D 26 are 
all zero within each plate strip, with the consequence that 


U R = V { = Wj = 0 


(B 5) 


so equation {B4) becomes 
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U(Yh 

r'U/OO') 


v(y) 

• = ^ R (V)l 

(B6) 

W(Y)J 

Iiv r (Y)J 



and the eigenfunctions (B3) can now be written in the form 
ru(x, yK r U,(Y) cos X 

J Kx,y) = V R (Y)sinXl (B7) 

U(x,y)J lkV R (Y)sinxJ 

In the following sections, the functions U,(Y), V*(Y), and W R (Y) are determined using the dif- 
ferential equations governing buckling. Symbols U, V, and W will be used henceforth to denote 
the functions U u V*. and W R , respectively. 


B.1.2 Out-of-Plane Displacements 


The out-of-plane portion of the buckling eigensolution 
ferential equation, taken from the body of the text: 

^X/'xx + 2M X y t ,xy + Myyyy ^i(^x L w l'xx ^y L w, yy) = ® 


is governed by the following dif- 


(B8) 


where unit loads N tL and N yL are known constants on a given plate strip, and the subscript /, 
which denotes the eigensolution index number, will henceforth be omitted. As shown in [5], 
through use of the assumed form for w(x, y) {equation (B7) ) and the definitions for the strains 
and stress resultants, the partial differential equation (B8) is converted to the following linear, 
homogeneous ordinary differential equation: 

\N"" - 2T1V" + ( T 2 - L)W = 0 (B9) 


where primes denote differentiation with respect to Y, and parameters T and L are given by 
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n 2 



T = ■ 


Ol2 + + A 


'22 


66 


^ 2 2 
20? rr 


L 2 N y 


L = ■ 


- X 


'22 


2 2 
m it 


± + t 2 -d 


11 


(BIO) 


The general solution for equation (B9) has the form 
W = Ce p/ (All) 

where C in an arbitrary constant, and p is a constant governed by the characteristic equation 
p 4 — 2Tp 2 + (T 2 — L) = 0 (512) 


or 


P 2 = T±7T (513) 

The exact nature of the roots p depends on the relative amplitudes of the parameters T and 
L, but the following common form will be used to express the function W(Y): 


W(Y)= (614) 

1 

where the constants C y are determined by the boundary conditions at the side edges of the 
plate strip. In [5], four cases are identified for characterizing the form of the function W(y), but 
the first case corresponds to conditions which are not considered here. The remaining three 
cases (Case (b) through Case (d)) are considered independently below. 

VIPASA Out-of*Plane Case b) L >0: Two non-negative, real roots are identified, consist- 
ent with equations (B13): 

a 2 = — oi 2 = T + y/U y 2 = -? 2 = r-7T (515) 

where either o or a is real, and either y ory is real. Note that if y is real then a must be real, 
and that if a is real then y must be real. Four roots p are then given by p = ± a or p = ± /a, 
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and p = ± y or p = ± iy. The four roots p establish four contributions to W{Y) t each having the 
form given in equation (B11), where it is understood that only the real portions of complex 
contributions are retained. The four contributions are combined and regrouped in such a way 
that the four functions fj of equation {B14) take the forms given in the following table; 

Table 6. Functional Forms Used for VIPASA Out-of-Plane Case (b) 



o 

Al 

« 2 >0 



o 

Al 

CM 

M 

V 

o 

'm 

COSh(aY) 

cos («Y) 


kY) 

cosh(yY) 

COS(y Y) 

h(Y) 

sinh(aY) 

sin (SlY) 


kY) 

sinh(yY) 

sin(yY) 


VIPASA Out-of-Plane Case c) L < 0: Here the roots of equation (B13) are complex, and 
two positive real constants a and /? are sought which satisfy 

p 2 = (« + /0) 2 = r±//T (si6) 

The values a and /? are found to satisfy 

« 2 = y(7+7 r2 - L ) >0 ^ 2 = y(-T + > /r 2 -L)>0 (SI 7) 


Four roots p are then given by p = ± (a ± //?). The four roots p establish four contributions to 
W(Y), each of the form given In equation (B11), where it is understood that only the real 
portions of complex contributions are retained. The four contributions are combined and re- 
grouped in such a way that the four functions fj of equation (B14) are given in the following 
table: 

Table 7. Functional Forms Used for VIPASA Out-of-Plane Case (c) 


U 

Y) 

cosh(aY) cos (/?Y) 

m 

Y) 

sinh(aY) cos(/?Y) 


Y) 

cosh(aY) sin (j?Y) 

m 

Y) 

sinh(aY) sin (/?Y) 



VIPASA Out-of-Plane Case d) L = 0: A single non-negative real value p or p is sought, 
defined by the equation 
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There are only two distinct roots, pi = p and Pi = — p, or pi = ip and p 2 = — ip, so the form of 
the solution for VV(V) is given by 

W(Y) = C 1 e p,v + C 2 e _PlV + C 3 ye PlV + C 4 ye _Piy (S19) 

where only the real portions of complex contributions are retained. The four contributions are 
combined and regrouped in such a way that the four functions fj of equation (B14) take the 
forms given in the following table: 

Table 8. Functional Forms Used for VIPASA Out-of-Plane Case (d) 



p 2 > 0 

£ 2 >0 

kY) 

cosh (pY) 

cos(£y) 

kY) 

sinh^y) 

y cos(£y> 

kY) 

ycosh^y) 

s\T)(pY) 

kY) 

y sinh(pV) 

Y sin(£y) 


B.1.3 In-Plane Displacements 


The in-plane portion of the buckling eigensolution is governed by the following equations, 
taken from the body of the text: 


^Xj’X + tyty/y ~ 0 

^xy/x T~ ^y,’y ^^x L v bxx ~ ® 


(B20) 


where the third term originally present in the first of these two equations has been omitted, 
consistent with the discussion in the main text, and the unit load N xL is known based on the 
prebuckling analysis. The subscript /, which denotes the eigensolution index number, will 
henceforth be omitted. Through use of the assumed forms for u(x,y) and v(x,y) of equations 
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(B7) and the definitions for the stress resultants, the partial differential equations (B20) can 
be converted to the following two coupled, linear, homogeneous, ordinary differential 
equations [5]: 


A 66 U" - L,U + A 0 V' = o 

—A 0 W + A 22 V" -L 3 V = 0 

where the newly introduced constants are defined by 

A 0 = ^12 + ^66 

l-l = A ii 

1 - 3 = ^66 + AW Xt 


(S21) 


(B22) 


Equations (B21) can be manipulated through differentiation and elimination of variables, with 
the result that both U and V must satisfy the same fourth-order differential equation: 


U"" - 2 BU" + CU = 0 
V" _ 2BV" + CV = 0 


(623) 


where constants 6 and C are given by 


6 = ■ 


^66^-3 "h AoyLj — A 


22H 
2A 22 A 66 


C = 


L<L 


1*-3 


^ 22^66 


(B 24) 


While U(Y) and V(Y) must individually satisfy their respective governing equations (B23), there 
is, in addition, a compatibility condition relating the two functions which can be determined 
using either of the equations (B21). 

The general form of the solutions for U and V are given by 

U(Y) = Ee pY V(Y) = Ge py (825) 


where E and G are arbitrary constants, and p is a constant governed by the characteristic 
equation 
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(B 26) 


p i -2Bp Z + C = 0 


or 

p 2 = e ± Ve 2 -c (627) 

The exact nature of the roots p depends on the relative amplitudes of the parameters B and 
C, but the following common form will be used to express the functions U{Y) and V(Y): 


U(V)=^EfV) 

/= 1 
4 

VM-^GjST/Y) 

/-I 

The compatibility condition can be expressed as 


(S28) 


£;= Yj R )* G * > =1 ' 2 ' 3 - 4 ( S29 ) 

k = 1 

where the values fy* are constants. The constants E t and G,of equations (B28) are determined 
by the boundary conditions at the side edges of the plate strip. In [5], two cases are identified 
for characterizing the form of the four functions g^Y). These cases are considered independ- 
ently below. 

VIPASA In-Plane Case a) B* > C: Two non-negative, real roots are identified, defined by 
the following equations: 

e 2 = - e 2 = b + ^/b 2 -c cf> 2 = -i> 2 = B- sJb 2 -c (B30) 

a A 

where either 9 or 6 is real, and either <f> or </> is real. Note that if <f> is real then 6 is real, and 

A A A 

similarly, if 9 is real, then <f> is real. Four roots p are given by p = ± 9 or p = ± id, and 
p = ± (f> or p = ± /(£. The four roots p establish four contributions to each of the functions 
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U(Y) and V(Y), consistent with the forms given in equations (B25), where it is understood that 
only the real portions of complex contributions are retained. The four contributions to each 
function are combined and regrouped into expressions having the form of equations (B28), and 
the compatibility constants R jk of equation (B29) are determined by substituting equations 
(B28) into either of the differential equations (B21). The four functions g) and the non-zero 
constants R Jk are given in the following table: 

Table 9. Functional Forms Used for VIPASA In-Plane Case (a) 



0 2 ^cr 

a, 

f? 2 > 0 



4> 7 ;>o* 

ro 

V 

0 

9 i(V) 

cosh(0V) 

cos(0y) 


9z(Y) 

cosh(^>y) 

cos{4>Y) 

9 2 (Y) 

sinh(0Y) 

sin(0V) 


94(Y) 

sinh(<£Y) 

s\n(4>Y) 

R 12 

R 

R 


R 34 

S 

s 

^21 

R 

- R 


^43 

s 

- $ 

* If e = 0: E, = E 2 = G, = G 2 = 0 
If $ * 0: Ej = £ 4 = G 3 = G 4 s 0 


The constants /?,R,S, and S introduced in the table above are given by 


R 


A 0 9 A 2 2& 2 ~~ I 3 

a 66 o 2 -l 1 " A ° 6 


R = 


f 

A 


A 0 e 


AqqO 2 + 


4- L 3 

A 

A o 0 


s = -- 


a 0 4> 

2 — 


A 7 2 < i > — 13 
A 0 <f> 


S = 


A 0 4> 

Agg4> 2 + L 1 


A22<f > + ^3 


A 0 4> 


(B 31) 


(B32) 


The two expressions given for each parameter above arise through separate applications of 
the first and second equilibrium equations (B21); the alternate forms can be shown to be 
equivalent. 

VIPASA In-Piane Case b) B 7 <> C: Here the roots p of equations (B27) are complex, and 
two non-negative real constants, <£1 and <f> 2 , are sought which satisfy 

p 2 = fa + i<j> 2 ) 2 = B± ijc - B 2 (S33) 
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The values <£i and <£ 2 are found to satisfy the equations 


<*M=y(V^ + B)>0 4>l = ±(jc -B)>0 (B34) 

Four roots p are then given by p = + ($1 ± /<£ 2 ). The four roots p establish four contributions 
to each of the functions U(Y) and Y(Y), consistent with the forms given in equations (B25), 
where it is understood that only the real portions of complex contributions are retained. The 
four contributions to each function are combined and regrouped into expressions consistent 
with equations (B28), and the compatibility constants R Jk are determined in the same manner 
as for In-plane Case (a). The four functions gj and the non-zero constants R Jk are given in the 
following table: 

Table 10. Functional Forms Used for VIPASA In-Plane Case (b) 



B 2 <C 

B=lVc'l (<f >2 = 0 ) 

s = — 1 VC 1 «M = 0 ) 

9i00 


cosh^Y) 

cos(4> 2 Y) 

9 2 (Y) 

sinh(<^ 1 V) cos(<£ 2 Y) 

sinh(<£.,y) 

Y cos(<j> 2 Y) 

9 3 (Y ) 

cosh (<f>^Y) sin(<£ 2 Y) 

Ycosh^Y) 

s\n{4> 2 Y) 

9a(Y) 

sinh(<^ 1 Y) sin {<f> 2 Y) 

Y sinh(<^ 1 Y) 

Y s\n(<f> 2 Y) 

P 12 . ^34 

P 

P 

P 


P 

P 

0 

1P1 

0 

Q 

Q 

Mjmxm 

-Q 

0 

-Q 


The constants P, P, Q, and Q are given by the following expressions: 


_ A 0 <t>i(A eey fc~ — Lj) <fti(Ag 2 \/C~ — L 3 ) 

(4 6 2 6 C - 2 A 66 L,B + L?) ~ A 0 jc 

— _ ^o(^66 ~ ^l) _ ~ ^3 

(A 66 Jc+L ,) 2 a oJc 
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_ 2(^66 V^~ + £1) ^2(^22%/^ + ^3) 

“ (4 6 2 6 C - 2A W L 1 B + L?) “ A 0 Jc 

— _ ^o(^ 66 %/^~ + £l) _ ^22 V^~ + ^-3 

(^ 66 Vc - Z.^ 2 


(836) 


The two expressions given for each parameter above arise through separate applications of 
the first and second equilibrium equations (B21); the alternate forms can be shown to be 
equivalent. 


B.1.4 Final Expression of the Eigensolutions 


It is desired to express the buckling eigensolutions in the form 

ru(x,y)-| f £(y) cos(mrtX/L)-| 

< V(x, y) > = < rj(y) s\rt(mnxlL) V (837) 

lw(x, y)) l<£(y) sin(m7tx/L)J 

where the subscript denoting the eigensolution index number is omitted, and tj, and <£ are 
given by 

(W) A fW] 

fiw =1 W ( S38 ) 

U(y)3 i-dc/jiy)) 

The barred and unbarred functions in equations (B14), (B28), and (B38) have the following re- 
lationship: 


U00] U(V)] U{mnylL)l 

(9y(y)j |g ; {V)j \gj(mnylL) J 


y = 1, 2. 3. 4 


(839) 
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Thus, in using the expressions for f t {Y) and gy(Y) presented in the appendix to express deriv- 
atives of the functions £,(y), r\{y), and <f>,(y), an appropriate constant must be incorporated, as 
shown in this example: 


di(y) y d W}_ mn y r 

dy ~ dy L /-j ‘ dY 

j = i ;= i 


( 840 ) 
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Appendix C. Finite-Difference Solution for the 
Second-Order Displacement Fields 


In this appendix, a procedure is established for determining an approximate finite- 
difference solution for the functions which privide the y-dependence of the second-order dis- 
placement fields. These functions are denoted as {£.//y)} = iLijiy) *\*ij{y) <^«->(y)] r . «= 1, 2 , 

and i,j = 1, 2 The solution procedure is performed once for each combination of a, /, and 

j (except for case ji when case ij has already been solved). With this understood, the sub- 
scripts a, /', and j will be omitted for the remainder of this discussion, except where needed for 
clarity. First, finite-difference approximations for derivatives are used to determine finite- 
difference expressions of the governing differential equations and the side-edge stress re- 
sultant contributions on a plate strip. A linear system of equations is formed in terms of the 
unknowns in the discretized y-domain of the plate strip, and then the appropriate transf- 
ormations are introduced which allow the assembly of a system of linear equations applying 
to the complete configuration. Boundary conditions are applied to the global system of 
equations. 


C.1.1 Finite-Difference Expressions for the Governing Equations 


On each plate strip, the functions {£} are governed by differential equations of the fol 
lowing form; 
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Cir + C 2 t + C&’ = F(y) 

D,? + D 2 >j" + *V? = G 00 

Erf’"' + E 2 4>" + E z 4> = H(y) 


(Cl) 


The associated generalized edge force-resultant amplitudes, used in expressing boundary 
conditions, depend on the functions {£} as shown in the following equations: 

^ = (-D e CCir + C 2 ^ + F(y)]|^ 

fy =(-V e lD,Z + D 2 r,' + G(y)]\ 

y ' e = 1, 2 (C2) 

4 =(-i) e [£ 1 ^"+E 2 ^ + «(y)^L 

e J9 

m e = -( -1) e [C 3 4>" + 1* <*3| 

y e 

The y-domain of the plate strip is discretized into / intervals, defining / + 1 discrete 
y-stations. The functional values at these discrete stations will be denoted by, for example, 
U 6, ... , {/. fi+i. Where 

tj-Wj) ( C3 > 

where the station ft is defined by 

y t -U- fl-f < c4 > 


where b is the width of the plate strip. 

Derivatives of the functions {£} at the discrete stations are expressed using finite- 
difference approximations. The formulae used are central-difference representations with 
truncation errors of order ( h 2 ) t where h is the interval width, given by 
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The finite-difference formulae are [63]: 


'■Ml, + + 

' n 

f "(y)\y=^-(- f j-2 + 2 f j-i-Vj + :+ f i + 2) 

r "0')ly=TT( / /-2-4/J-i + 6f y -4f ; + 1 + f 7 + 2 ) 

y h 


Using the above formulae for derivatives to express the differential equations (Cl) at y = y y , 
the following finite difference equations are obtained: 

^1l£y- 1 + ^12 £j + ^13^ + 1 + k u y lj- 1 + +1 — Fj 

k 21^j — 1 + k 22^j + 1 + k 23*lj — 1 + k 24Vj + k 25*1j + 1 = J = 2 . 1,1+ ^ (C7) 

*31^) -2 + ^32^; -1 + ^33^y + ^34$y + 1 + ^35^y + 2 = 

where the constant coefficients are given by 



and the right-hand-side functions are given by 

F y = F(y y ) G y = G(y y ) « y = H(y y ) (C9) 

It is helpful to express equations (C7) in terms of matrices and vectors. Introduce the 
notation {{ y } = [{>»?/ <f>jT and {F y } = [F y G y H y ] T , and now equations (C7) can be written as 

[K- 2 ]{£ y _ 2 } + [K -1 ]{{y_ 1 } + [K°]{{ y } + [K +1 ]{{ y + 1 } + [K +2 ]{(* y + 2 } = {F y } (CIO) 
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where 


* 13*15 0 
*22 *25 0 
0 0 

and where the only nonzero elements of [K -2 ] and [K* 2 ] are given by 
K33 = *31 

1/+2 _ i, 
a 33 “ k 35 


’*11 

*14 0 


1 

_ » 
ro 

O 

O 


/( 2 1 

*23 0 

[K°] = 

0 k 24 0 

II 

1 — 1 

+ 

* 
1 1 

0 

0 A 32 j 


0 0 ^33 



(C11) 


(C12) 


C.1.2 Finite-Difference Expressions for the Generalized Edge 
Force-Resultants 


As with the governing equations, the edge stress-resultant amplitudes are expressed in 
terms of discrete values of the displacement variables by applying the finite-difference for- 
mulae of equations (C5) to equations (C2). The following expressions are obtained: 


£., = 4^0 -|^ 2 -C 2 >h-F(0) 

• — Do Do — 

k, = - £1 + >fo - ^ ^2 - G (°) 




<^3-H(0) 


(Cl 3) 
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Cl *' + |h' + 2 + C 2 f, l + , + F(b) 


2h 


f y ' 2 = D^ l + ,- — > 1 i + — >ii + 2 + G(b) 


f i*2~ _ . ^ — i + ( i3 2 h i ^ 1 + » . 3 + ?h r 1 


m„ = — ■ 


2h' 


2h 


1 + 2 


2h' 


■4>, + 3 + H(b) 


(C14) 


4>i + 


2£ 


3 c 


-e 4 


+ 2 


where ei implies e = 1, and e 2 implies e = 2. 

Inspection of equations (C7), (C13), and (C14) reveals that eight nonexistent function val- 
ues have been introduced, corresponding to stations outside of the y-domain. These are de- 
noted as <p^, U r]o, and <f> a , associated with values of y less than zero, and 
{/ + *. n '»+*, 4*i + 2 i and <f> l + 3 , associated with values of y greater than b. It is necessary to replace 
these nonexistent values with expressions involving existent function values. The first step in 
doing this is to add an edge-rotation variable to each end of the y-domain, denoted and 
defined by 


k = 4>'(y) I. 


e = 1,2 


(Cl 5) 


This establishes 3(/+ 1) + 2 discrete unknowns across the width of the plate strip. Recalling 
that {£,} = [£, i], 4> t fa] (e = 1, 2) and introducing the notation 

{cf} = [£ 2 m 4 > 2 & »/3 4*3 ■■■ Zt Vi 4>.T, the vector of unknowns is given by 


r{t*h 

{d} 

^ et y 


(C 16) 


Express in terms of finite-difference approximations to get 




xt/ ^ = ~~2h 4, ‘ + ~2h 4 ‘ l + 2 


(C 17) 
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By rearranging equations (C17), the following expressions for <f> 0 and <£ / + 2 can be obtained: 


4> 0 = — + 4> 2 <£/ + 2 “ ^ x l / e 2 + <£/ (Cl®) 

Equations (C18) are also written another way for later use: 

© 0 

<^0 = C 1^1 +C 2 4>2 $1 + 2 ~ c 1^2 + c 2^/ (C19) 

where the constant coefficients are given by 
c 1 = — 2h c 2 = 1 

(C20) 

c^ = 2h c 2 = 1 

The six remaining nonexistent function values can be expressed in terms of existent val- 
ues by rearranging equations (C13) and (C14), and making substitutions using equations (C18) 
where appropriate. The following expressions can be obtained: 

£/ + 2 = 3/xe* + ^2 £/ + + 1 + 

*l/ + 2 = 1 + ^3*?/ + ^4 (C21) 

_ • __ © 

<I>I + 3 = d^f„ 2 + d 2 4 ! et + d z 4> I _ , + d 4 

where 


(C 22) 


Appendix C. Finite-Difference Solution for the Second-Order Displacement Fields 194 


t> 2 = 2h^~ 
D 2 


2 h 3 E 2 


d 0 = -4h + — =7 


a z = 2h^ 

c i 


b 3 = 1 


c/, = 1 


3 H(0) 
cf 4 = 2/T — =-^- 


0 

fo = a 1^ei + a 2^2 + a 3*?1 + a 4 
« 

>fo = + ^2^1 + ^3^2 + ^4 

<f >- 1 = + d z <f> z + d 4 



(C 23) 


2h 

Ci 

a 2 — 1 

c 2 

a 3 = — 2h 

c i 

a 4 — 

F(b) 

- 2h 

Ci 

2h_ 

d 2 

b 2 = - 2h -5- 

b 3 =1 

b 4 = 

G(b) 

-2h^ J - 

d 2 

2h 3 

r~ 

2 h 3 E 2 

d 7 = Ah- _ 

* r 

d 3 = 1 

i 5 -' 

II 

3 H(b) 

— 2h 3 -=r$- 
c 


The edge moments of equations (C13) and (C14) were not required in determining the 
expressions of equations (C21), but they will be required in order to add two additional 
equations corresponding to the two edge-rotation amplitudes which have been added as un- 
knowns. For this purpose, the edge moments of equations (C13) and (C14) are expressed in 
the following way: 


^ , e 1 = 01^0, + 92 $ 1 + 93$ 2 


m e 2 ~ 9'i$e z + 92$ l + 93$ I + i 


(C24) 


where 


01 = - 





01 = - 


2E3 

h 




(C 25) 


Some new notation Is introduced in order to facilitate expression of equations (C19) and 
{C21) in terms of vectors and matrices. Define the vectors {{»} and {£*} (e = 1, 2): 



(C26) 


Expressing {£,} (e = 1, 2) in terms of finite-difference unknowns and \j/„ we can say that 
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(C 27) 




£1 

£/ + 1 

s 

fn + i 
Ke 2 } = 1 

4> ij 

U/+1 

O 

0 


'l'e 2 


Now equations (Cl 9) and (C21) can be expressed using matrix notation: 


{Q = !> 0 ]{V} + O 1 ^,} + l> 2 ]K 2 } + {^ 4 } 

= [S°]{^/} + [B 1 ]^,} + [S 3 ]R 3 } + {B 4 } 


{£, +2 } = C^ 0 ](V} + + (* 4 } 

{e /+3 } = Cb°]{^*} + [ri{i e2 } + Cs 3 ]{^_ 1 } + {b 4 } 


(C 28) 


(C29) 


The nonzero elements of the matrices introduced In equations (C28) are given by 


a 0 i 0 i. 

^11 = a 1 ^22 “ 

^12 = a 3 ^21 = ^2 ^34 ~ C 1 

2 2 2 
^11 = a 2 ^22 = ^3 ^33 = ^2 

/A >| = a 4 — b 4 

The nonzero elements of the matrices in equations (C29) are given by the same expressions 
as in the equations above, except that the unsubscripted constants a, b, and d are replaced 
by the subscripted constants a, b , and d. Now ail the elements are in place for assembly of 
the system of finite-difference equations which apply to a single plate strip. 


e 33 = d 1 


^34 “ d 2 


S33 — d 3 

B 3 4 = d 4 


(C30) 
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C.1.3 System of Equations for a Plate Strip 

Equations at y=y,: To obtain the first three equations, consider the finite-difference 
equations (CIO) for the case j= 1 {y - 0): 

!*-*]{«_,} + [K -1 ]{{<,} + [K°]{| ei } + [K +1 ]{^} + [K +2 ]{£ 3 } = (F-i) (C31) 

where {£«,} of equation (C26) has replaced {£ 1 }, and one coefficient matrix was modified to 
account for the added degree of freedom: 

I o' 

[K°] = [K°] | 0 (C32) 

I 0 

The vectors {<*;_,} and {£<,} appearing in equations (C31) are replaced by the expressions given 
in equations (C28); equations (C31) can then be expressed in the form 

[*]&,} + u]{£ 2 } + [«]{« = Mv) + {^i> ( C33 ) 

where the newly introduced matrices and vector are given symbolically, and in evaluated 
form, by 

(/C12 + ^14^2) ^11 a 3 0 0 

£/(] = [K°] 4 CK ] 4 2 ][e ] = ^23^2 (^24 + ^21 a 3) ® ® (C34) 

0 0 k 33 (/f32Ci 4 ^ 31 ^ 2 ) 

(^13 + ^11 a 2)(^15 + *14^3) 0 

[/- ] = [K +1 ] 4- [K 1 D[^A = (^22 + ^2l a 2)(^25 + ^ 23 ^ 3 ) ® (C35) 

0 0 (JC 34 4 ^ 32 c 2 ) 
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[J] = - ([k - 1 ][/4°] + [K~ z ][e 0 ]) = - 


*21 a 1 * 23^1 ^ 

0 0 *31 


Fi — (^ii a 4 + *1464)^ 


{£,} = (F,} - ([K -1 ]^ 4 } + [K- z ]{8 4 }) = {G,~ (k 21 a 4 + k 23 b A ) 


— 2 - 


H-i ~ k^d 4 


and where [/W] has only one nonzero element: 

C/W] = [K +2 ] + [K 2 ][S 3 ] , — (^35 + ^ 31 ^ 3 ) 

Finally, isolate {4,*} by premultiplying equation (C33) by [J] 1 : 

[K ]{?.,} + CH{« + [M]{«3) = {4;} + {Fi} 


where 


([K], [L], [m3. {F}) = [J] \ck], [l], [m3. (fj) 


(C36) 


(C37) 


(C38) 


(C39) 


(C40) 


The three linear equations (C39) correspond to the first three unknowns, \, v f) ev and of the 
vector in equation (C16). 

The fourth equation, corresponding to the unknown ^ ei , is the first of equations (C24). 
Rewriting the equation in notation compatible with equations (C39), 


[0 0 6^2 + [00 s ? 33 {<^2} — 


(C41) 


Equations at y = y 3 : To obtain the fifth through seventh equations, corresponding to un- 
knowns It, >j 2 , and 4 > 2 of the vector in equation (C16), consider the finite-difference equations 
(CIO) for the case j - 2: 
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[K -2 ]{£ 0 } 4- iK-'lile,) + [X°]{^ 2 } + [K +1 ]{£3} + [K +2 ]{^ 4 } = {^} 


(C 42) 


where {£„} of equation (C26) has replaced {£ 1 }, and one coefficient matrix was modified to 
accommodate the added degree of freedom: 


DC 1 ] 


[K _1 ] I 0 
I o 


(c 43) 


The vector {£„} appearing in equations (C42) is replaced by the expression given in equations 
(C28). Equations (C42) can then be expressed in the following form, which is suitable for direct 
incorporation into the system of equations for the plate strip as a whole: 


+ [k°]{£ 2 } + [k +1 ]{£ 3 ) + Ck +2 ]{£ 4 } = { f 2 } 


(C 44) 


where the newly introduced matrices and vector are given symbolically, and in evaluated 
form, by 


[k- 1 ] = [k- 1 ] + [k- 2 ][/i 1 ] = 


-2-ir*i- 


*11 *14 0 


0 

0 0 
0 0 /c 32 ^3-|C-| 


*21 *23 


(C 45) 


[K°] = [K°] + [K- 2 ][/t 2 ] = 


0 0 

0 ^24 0 

0 0 (fc 3 3 + ^31 C 2 ) 


(C46) 


and where it is noted in addition that some expected matrix terms do not appear in equations 
(C44) because certain matrix products vanish: 

[K- 2 ][4°] = { 0} [k- 2 ]{4 4 } = { 0} (C47) 
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Equations at y= y,: The three equations evaluated at y, require the same special con- 
sideration as those evaluated at y z . Consider the finite-difference equations (CIO) for the case 
J = /: 


[K- 2 ]{£,_ 2 } + [/C 1 ]^} + [K° ]{£,} + [K +1 ]{^} + [K +2 ]{{, + 2 } = {F 2 } (C48) 


where of equation (C26) has replaced {{ / + 1 }, and one coefficient matrix was modified to 
accommodate the added degree of freedom: 


[K +1 ] = 


10 

[K +1 ] I 0 

I 0 


(C49) 


The vector {& + ?} appearing in equations (C48) is replaced by the expression given in 
equations (C29). Equations (C48) can then be expressed in the following form, which is suit- 
able for direct incorporation into the system of equations for the plate strip as a whole: 


[K- 2 ]{£,_ 2 } + [ + [K°]{£,} + [K +1 ]{^} = {F,} 


(C50) 


where 


[K°] = [K°] + [K +2 ][^] = 


/ c<!2 0 0 

0 ^24 0 

0 0 (ff 3 3 + ((35 C 2 ) 


(C51) 


[K +1 ] = [K +1 ] + [k + 2 ][V] = 


*13 *15 0 0 

*22 *25 ^ 0 
0 0 /C34 JC35C.1 


(C 52) 


and where it Is noted in addition that some expected matrix terms do not appear in equations 
(C50) because certain matrix products vanish: 

[K +2 ][A°] = {0} [K +2 ]{* 4 } = { 0} (C53) 
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Equations at y = y, + 1 : The three equilibrium equations evaluated aty / + 1 require the same 
special consideration as those evaluated at yi. Consider the finite-difference equations (CIO) 
for the case j = I + 1 (y = b): 

[K" 2 ]K,_. l) + [K" 1 ]^} + CK 0 ]{| e2 } + [K +1 + + [K +2 ]{{ /+3 } = (F, + 1 ) (C54) 

where {£ ej } of equation (C26) has replaced {£ /+1 }, and the modified coefficient matrix [K°] is 
given by equation (C33). The vectors {{/ + *} and {& + 3 } appearing in equations (C54) are re- 
placed by the expressions given in equations (C30); equations (C54) can then be expressed 
in the form 

[#]{*/_ i) + [£]{*,} + [K]{^} = [J]{4/} + £>} (C55) 

where the newly introduced matrices and vector are given symbolically, and in evaluated 
form, by 


[K] = [K°] + [K +1 ][^] + [K +2 ] tr ] = 


(*12 + * 15 ^ 2 ) * 13 a 3 ^ 0 

* 25^2 (*24 + * 22 a 3 ) ^ ® 

0 0 *33 (*340-1 + *35^2) 


(C56) 


[L] = [K- 1 ] + [K +1 ][^] = 


(*11 + *13 a 2)(*14 + *15^3) 0 

(*21 + *22 a 2)(*23 + *25^3) 0 

0 0 (*32 + *34^2) 


(C57) 


* 13 a 1 * 15^1 0 


[J] = - ([K +1 ][^°] + [K +2 ][B 0 ]) = - 


*22 a 1 *25^1 0 


0 0 * 35^1 


(C 58) 
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/) + 1 “ ( fc 13 a 4 + *15 b 4). 


{F / + 1 } = {F, + 1 }- ([K 1 ]^ 4 } + [K 2 ]{6 4 }) = <| G, + 1 - (* 22 a 4 + * 25 b 4 ) 


(C59) 


Hi + 1 _ * 35^4 


and where [/W] has only one nonzero element: 

[M] = [K ^3 + [K + 2 ][S 3 D , M 33 = (/C 3 1 + /C 35 C/ 3 ) (C60) 

Finally, isolate {£/} by premultiplying equation (C55) by [J ] -1 : 

Cm]«,_i} + [£]{«,} + = {^‘} + {F, + 1 } (C61) 

where 

([K], [Q [M], {F}) = [J ] _1 ([i<], [L], [M], {F}) (C62) 

• « 

The three linear equations (C61) correspond to the unknowns r\ 9V and 4> ez 

The equation corresponding to the unknown \j/ e2 is the second of equations (C24). Re- 
writing the equation in notation compatible with equations (C61), 

[0 0 + tO 0 03 = m e2 (C63) 

Equations at y = y 3l y 4 , ... For the y-stations corresponding to j = 3, 4 / - 1 , the 

application of equation {CIO) in assembling the stiffness matrix for the plate strip is 
straightforeward, with the one clarification that for/’ =3 andy = /— 1 , the vectors {{ 1 } and 
{£, + 1 } which appear are expressed in terms of the vectors {£ ei } and {£ ez }, respectively (see 
equation (C27)). 

Note regarding the preceding derivations: Close Inspection of the development in this 
section, beginning with equation (C31), will reveal that the matrix operations as used were 
unnecessarily cumbersome, because the third of the governing equations (C7) remains un- 
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coupled from the first two throughout the derivations. Nonetheless, the development is correct 
as presented. 

Final expression of the equations for a plate strip: Now a system of linear equations has 
been assembled which applies to the plate strip. Introducing the notation 
[F] = [f 2 G 2 H z F 3 G 3 H 3 ... F, Gi H/Y, and recalling the notation for the displacement unknowns 
given in equation (C16), the system of equations can be expressed in the following way: 


*11 I K12 I o {l e J (f e ) (F,) 
K 2 i I K 2 j I K 23 {cf} = {0} + { F } 


0 I K 32 I K 33 


«e*> {4 2 } {F 2 } 


(C64) 


where certain sub-matrices are zero as indicated, assuming / is greater than or equal to four. 

In order to assemble the systems of equations for the individual plate strips into a global 
system of equations for the entire panel, the vectors for generalized edge displacements and 
generalized edge force-resultants, {£,} and {£}, respectively, must be transformed into the 
equivalent values referred to the appropriate node-line and the global coordinate axes. The 
transformations are effected through the following equations: 


tie) = lT ecc -] e lT r -]{U n '} 

e = 1, 2 

(C65) 

{f n *) = Cr r ] T [r ecc ]^; e } 

e = 1,2 

(C66) 


where the transformation matrices are given in the main text. To incorporate the transf- 
ormations into the system of equations, replace {£,} (e=1,2) with the appropriate ex- 
pressions given in equation (C65), premuitiply the first four equations of (C64) by 
and premultiply the last four equations of (C64) by [r,] T [T ecc ], 2 Now the system 
of equations has the form 
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0 


K 11 I K 12 


K21 I K< 


22 1 ^23 


^32 I ^33 


|{i) ni } (f^ {FO 

{of} ={0}+ {F} 

{u" 2 } {f" 2 ! {f 2 } 


where 

[k„] - [r r ] r [r ee{ .]^[K 11 ][r,„] e [7 r ] 
CK 12 ]-[7 r ]'[r,„] r [K, 2 ] 

[K 21 ] = [K 21 ][r ecc ] ei [r r ] 

[k 23 ] = [K 23 ][r ecc ]^[T,.] 

[K 32 ] = [F r ] T [r ecc ]^[K 32 ] 

[K 33 ] = Cr r ] T [7 ecc ]^[K3 3 ]Cr ecc ] ej [T r ] 
{F 1 } = [r f ] r [7 ecc ]^{F 1 } 

{F 2 } - [r r ] r [r ecc ]^{F 2 } 


(C67) 


(C68) 


C.1.4 Global System of Equations 


The global system of equations, applying to the entire panel, has as unknowns the four 
components of { U n } for each node line, and the 3 (/ - 1) values in { d) for each plate strip. When 
the system of equations applying to each plate strip (equations (C67)) has been assembled 
into the global system of equations, the sum of the vectors (f 1 ) for all plate edges terminating 
at a given node line is simply the vector of generalized force-resultant amplitudes for the node 
line, {A 1 }. Along nonboundary node lines, these values must be zero. Along boundary node 
lines, the homogeneous form of the selected options for boundary conditions are applied, so 
that each component of {A} either is zero or has its corresponding generalized displacement 
component set to zero, removing from consideration the equation containing the generalized 
force-resultant amplitude. Thus, the final system of equations to be solved after application 
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of boundary conditions has a form analogous to that of equation (C67), but without the pres- 
ence of the first right-hand-side vector. The stiffness matrix will in general be banded, and this 
fact is used to select an efficient numerical procedure for solution of the system of equations. 

Solution of the global system of equations provides the values values {()"} for each node 
line, and the 3(1 - 1) values in {d} for each plate strip. Using the transformation relationship 
of equation (C65), the vectors {£,} can be determined for each edge of each plate strip, and 
that completes the set of discrete values for (£(y)} on the plate strip. Also needed for later use 
are first and second derivatives of £(y) and tj (y), and first through fourth derivatives of <£(y). In 
order to determine these derivatives, the following sequence of operations is performed: 

1. Apply equations (C39) and (C61) to solve for vectors {£*} (e= 1, 2). 

2. Apply equations (C28) and (C29) to determine the nonexistent function values 
{f-i}. {«.{&♦«}. and {ft*,}. 

3. Use equations (C6) to determine finite-difference approximations for the derivatives at the 
end points of the y-domain. 
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